source: palm/trunk/SOURCE/microphysics_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

  • Property svn:keywords set to Id
File size: 86.6 KB
Line 
1!> @file microphysics_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21! Module renamed
22! Adapted for modularization of microphysics.
23!
24! Former revisions:
25! -----------------
26! $Id: microphysics_mod.f90 1850 2016-04-08 13:29:27Z maronga $
27!
28! 1845 2016-04-08 08:29:13Z raasch
29! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
30! point (instead of one point above surface)
31!
32! 1831 2016-04-07 13:15:51Z hoffmann
33! turbulence renamed collision_turbulence,
34! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
35! avaialble for microphysics_kessler.
36!
37! 1822 2016-04-07 07:49:42Z hoffmann
38! Unused variables removed.
39! Kessler scheme integrated.
40!
41! 1691 2015-10-26 16:17:44Z maronga
42! Added new routine calc_precipitation_amount. The routine now allows to account
43! for precipitation due to sedimenation of cloud (fog) droplets
44!
45! 1682 2015-10-07 23:56:08Z knoop
46! Code annotations made doxygen readable
47!
48! 1646 2015-09-02 16:00:10Z hoffmann
49! Bugfix: Wrong computation of d_mean.
50!
51! 1361 2014-04-16 15:17:48Z hoffmann
52! Bugfix in sedimentation_rain: Index corrected.
53! Vectorized version of adjust_cloud added.
54! Little reformatting of the code.
55!
56! 1353 2014-04-08 15:21:23Z heinze
57! REAL constants provided with KIND-attribute
58!
59! 1346 2014-03-27 13:18:20Z heinze
60! Bugfix: REAL constants provided with KIND-attribute especially in call of
61! intrinsic function like MAX, MIN, SIGN
62!
63! 1334 2014-03-25 12:21:40Z heinze
64! Bugfix: REAL constants provided with KIND-attribute
65!
66! 1322 2014-03-20 16:38:49Z raasch
67! REAL constants defined as wp-kind
68!
69! 1320 2014-03-20 08:40:49Z raasch
70! ONLY-attribute added to USE-statements,
71! kind-parameters added to all INTEGER and REAL declaration statements,
72! kinds are defined in new module kinds,
73! comment fields (!:) to be used for variable explanations added to
74! all variable declaration statements
75!
76! 1241 2013-10-30 11:36:58Z heinze
77! hyp and rho have to be calculated at each time step if data from external
78! file LSF_DATA are used
79!
80! 1115 2013-03-26 18:16:16Z hoffmann
81! microphyical tendencies are calculated in microphysics_control in an optimized
82! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
83!
84! 1106 2013-03-04 05:31:38Z raasch
85! small changes in code formatting
86!
87! 1092 2013-02-02 11:24:22Z raasch
88! unused variables removed
89! file put under GPL
90!
91! 1065 2012-11-22 17:42:36Z hoffmann
92! Sedimentation process implemented according to Stevens and Seifert (2008).
93! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
94! and Stevens, 2010).
95!
96! 1053 2012-11-13 17:11:03Z hoffmann
97! initial revision
98!
99! Description:
100! ------------
101!> Calculate bilk cloud microphysics.
102!------------------------------------------------------------------------------!
103 MODULE microphysics_mod
104
105    USE  kinds
106
107    IMPLICIT NONE
108
109    LOGICAL ::  cloud_water_sedimentation = .FALSE.  !< cloud water sedimentation
110    LOGICAL ::  limiter_sedimentation = .TRUE.       !< sedimentation limiter
111    LOGICAL ::  collision_turbulence = .FALSE.       !< turbulence effects
112    LOGICAL ::  ventilation_effect = .TRUE.          !< ventilation effect
113
114    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
115    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
116    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
117    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
118    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
119    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
120    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
121    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
122    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
123    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
124    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
125    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
126    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
127    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
128    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
129    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
130    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
131    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
132    REAL(wp) ::  eps_sb = 1.0E-20_wp       !< threshold in two-moments scheme
133    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
134    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
135    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
136    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
137    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
138    REAL(wp) ::  kappa_rr = 60.7_wp        !< const. in collision kernel (kg-1/3)
139    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
140    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
141    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
142    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
143    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
144    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
145    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
146    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
147    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
148    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
149
150    REAL(wp) ::  c_sedimentation = 2.0_wp  !< Courant number of sedimentation process
151    REAL(wp) ::  dpirho_l                  !< 6.0 / ( pi * rho_l )
152    REAL(wp) ::  dt_micro                  !< microphysics time step
153    REAL(wp) ::  nc_const = 70.0E6_wp      !< cloud droplet concentration
154    REAL(wp) ::  dt_precipitation = 100.0_wp !< timestep precipitation (s)
155    REAL(wp) ::  sed_qc_const              !< const. for sedimentation of cloud water
156    REAL(wp) ::  pirho_l                   !< pi * rho_l / 6.0;
157
158    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nc_1d  !<
159    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nr_1d  !<
160    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_1d  !<
161    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qc_1d  !<
162    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qr_1d  !<
163    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_1d   !<
164
165    SAVE
166
167    PRIVATE
168    PUBLIC microphysics_control, microphysics_init
169
170    PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation,   &
171           dt_precipitation, limiter_sedimentation, nc_const, sigma_gc,        &
172           ventilation_effect
173
174    INTERFACE microphysics_control
175       MODULE PROCEDURE microphysics_control
176       MODULE PROCEDURE microphysics_control_ij
177    END INTERFACE microphysics_control
178
179    INTERFACE adjust_cloud
180       MODULE PROCEDURE adjust_cloud
181       MODULE PROCEDURE adjust_cloud_ij
182    END INTERFACE adjust_cloud
183
184    INTERFACE autoconversion
185       MODULE PROCEDURE autoconversion
186       MODULE PROCEDURE autoconversion_ij
187    END INTERFACE autoconversion
188
189    INTERFACE autoconversion_kessler
190       MODULE PROCEDURE autoconversion_kessler
191       MODULE PROCEDURE autoconversion_kessler_ij
192    END INTERFACE autoconversion_kessler
193
194    INTERFACE accretion
195       MODULE PROCEDURE accretion
196       MODULE PROCEDURE accretion_ij
197    END INTERFACE accretion
198
199    INTERFACE selfcollection_breakup
200       MODULE PROCEDURE selfcollection_breakup
201       MODULE PROCEDURE selfcollection_breakup_ij
202    END INTERFACE selfcollection_breakup
203
204    INTERFACE evaporation_rain
205       MODULE PROCEDURE evaporation_rain
206       MODULE PROCEDURE evaporation_rain_ij
207    END INTERFACE evaporation_rain
208
209    INTERFACE sedimentation_cloud
210       MODULE PROCEDURE sedimentation_cloud
211       MODULE PROCEDURE sedimentation_cloud_ij
212    END INTERFACE sedimentation_cloud
213 
214    INTERFACE sedimentation_rain
215       MODULE PROCEDURE sedimentation_rain
216       MODULE PROCEDURE sedimentation_rain_ij
217    END INTERFACE sedimentation_rain
218
219    INTERFACE calc_precipitation_amount
220       MODULE PROCEDURE calc_precipitation_amount
221       MODULE PROCEDURE calc_precipitation_amount_ij
222    END INTERFACE calc_precipitation_amount
223
224 CONTAINS
225!------------------------------------------------------------------------------!
226! Description:
227! ------------
228!> Initialization of bulk microphysics
229!------------------------------------------------------------------------------!
230    SUBROUTINE microphysics_init
231
232       USE arrays_3d,                                                          &
233           ONLY:  dzu
234
235       USE constants,                                                          &
236           ONLY:  pi
237
238       USE cloud_parameters,                                                   &
239           ONLY:  rho_l
240
241       USE control_parameters,                                                 &
242           ONLY:  microphysics_seifert
243
244       USE indices,                                                            &
245           ONLY:  nzb, nzt
246
247       IMPLICIT NONE
248
249!
250!--    constant for the sedimentation of cloud water (2-moment cloud physics)
251       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
252                          )**( 2.0_wp / 3.0_wp ) *                             &
253                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
254
255!
256!--    Calculate timestep according to precipitation
257       IF ( microphysics_seifert )  THEN
258          dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
259                             w_precipitation
260       ENDIF
261
262!
263!--    Pre-calculate frequently calculated fractions of pi and rho_l
264       pirho_l  = pi * rho_l / 6.0_wp
265       dpirho_l = 1.0_wp / pirho_l
266
267!
268!--    Allocate 1D microphysics arrays
269       ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
270                  qc_1d(nzb:nzt+1) )
271
272       IF ( microphysics_seifert )  THEN
273          ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
274       ENDIF
275
276!
277!--    Initialze nc_1d with nc_const
278       nc_1d = nc_const
279
280    END SUBROUTINE microphysics_init
281
282
283!------------------------------------------------------------------------------!
284! Description:
285! ------------
286!> Control of microphysics for all grid points
287!------------------------------------------------------------------------------!
288    SUBROUTINE microphysics_control
289
290       USE arrays_3d,                                                          &
291           ONLY:  hyp, pt_init, prr, zu
292
293       USE cloud_parameters,                                                   &
294           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
295
296       USE control_parameters,                                                 &
297           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
298                  intermediate_timestep_count, large_scale_forcing,            &
299                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
300                  pt_surface, rho_surface,surface_pressure
301
302       USE indices,                                                            &
303           ONLY:  nzb, nzt
304
305       USE kinds
306
307       USE statistics,                                                         &
308           ONLY:  weight_pres
309
310       IMPLICIT NONE
311
312       INTEGER(iwp) ::  k                 !<
313
314       REAL(wp)     ::  t_surface         !<
315
316       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
317!
318!--       Calculate:
319!--       pt / t : ratio of potential and actual temperature (pt_d_t)
320!--       t / pt : ratio of actual and potential temperature (t_d_pt)
321!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
322          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
323          DO  k = nzb, nzt+1
324             hyp(k)    = surface_pressure * 100.0_wp * &
325                         ( ( t_surface - g / cp * zu(k) ) /                    &
326                         t_surface )**(1.0_wp / 0.286_wp)
327             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
328             t_d_pt(k) = 1.0_wp / pt_d_t(k)
329             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
330          ENDDO
331
332!
333!--       Compute reference density
334          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
335       ENDIF
336
337!
338!--    Compute length of time step
339       IF ( call_microphysics_at_all_substeps )  THEN
340          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
341       ELSE
342          dt_micro = dt_3d
343       ENDIF
344
345!
346!--    Reset precipitation rate
347       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
348
349!
350!--    Compute cloud physics
351       IF ( microphysics_kessler )  THEN
352
353          CALL autoconversion_kessler
354          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
355
356       ELSEIF ( microphysics_seifert )  THEN
357
358          CALL adjust_cloud
359          CALL autoconversion
360          CALL accretion
361          CALL selfcollection_breakup
362          CALL evaporation_rain
363          CALL sedimentation_rain
364          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
365
366       ENDIF
367
368       CALL calc_precipitation_amount
369
370    END SUBROUTINE microphysics_control
371
372!------------------------------------------------------------------------------!
373! Description:
374! ------------
375!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
376!> evaporation of rain drops due to too small or too big weights
377!> of rain drops (Stevens and Seifert, 2008).
378!------------------------------------------------------------------------------!
379    SUBROUTINE adjust_cloud
380
381       USE arrays_3d,                                                          &
382           ONLY:  qr, nr
383
384       USE cloud_parameters,                                                   &
385           ONLY:  hyrho
386
387       USE cpulog,                                                             &
388           ONLY:  cpu_log, log_point_s
389
390       USE indices,                                                            &
391           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
392
393       USE kinds
394
395       IMPLICIT NONE
396
397       INTEGER(iwp) ::  i                 !<
398       INTEGER(iwp) ::  j                 !<
399       INTEGER(iwp) ::  k                 !<
400
401       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
402
403       DO  i = nxl, nxr
404          DO  j = nys, nyn
405             DO  k = nzb_s_inner(j,i)+1, nzt
406                IF ( qr(k,j,i) <= eps_sb )  THEN
407                   qr(k,j,i) = 0.0_wp
408                   nr(k,j,i) = 0.0_wp
409                ELSE
410                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
411                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
412                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
413                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
414                   ENDIF
415                ENDIF
416             ENDDO
417          ENDDO
418       ENDDO
419
420       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
421
422    END SUBROUTINE adjust_cloud
423
424
425!------------------------------------------------------------------------------!
426! Description:
427! ------------
428!> Autoconversion rate (Seifert and Beheng, 2006).
429!------------------------------------------------------------------------------!
430    SUBROUTINE autoconversion
431
432       USE arrays_3d,                                                          &
433           ONLY:  diss, dzu, nr, qc, qr
434
435       USE cloud_parameters,                                                   &
436           ONLY:  hyrho
437
438       USE control_parameters,                                                 &
439           ONLY:  rho_surface
440
441       USE cpulog,                                                             &
442           ONLY:  cpu_log, log_point_s
443
444       USE grid_variables,                                                     &
445           ONLY:  dx, dy
446
447       USE indices,                                                            &
448           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
449
450       USE kinds
451
452       IMPLICIT NONE
453
454       INTEGER(iwp) ::  i                 !<
455       INTEGER(iwp) ::  j                 !<
456       INTEGER(iwp) ::  k                 !<
457
458       REAL(wp)     ::  alpha_cc          !<                   
459       REAL(wp)     ::  autocon           !<
460       REAL(wp)     ::  dissipation       !<
461       REAL(wp)     ::  k_au              !<
462       REAL(wp)     ::  l_mix             !<
463       REAL(wp)     ::  nu_c              !<
464       REAL(wp)     ::  phi_au            !<
465       REAL(wp)     ::  r_cc              !<
466       REAL(wp)     ::  rc                !<
467       REAL(wp)     ::  re_lambda         !<
468       REAL(wp)     ::  sigma_cc          !<
469       REAL(wp)     ::  tau_cloud         !<
470       REAL(wp)     ::  xc                !<
471
472       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
473
474       DO  i = nxl, nxr
475          DO  j = nys, nyn
476             DO  k = nzb_s_inner(j,i)+1, nzt
477
478                IF ( qc(k,j,i) > eps_sb )  THEN
479
480                   k_au = k_cc / ( 20.0_wp * x0 )
481!
482!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
483!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
484                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
485!
486!--                Universal function for autoconversion process
487!--                (Seifert and Beheng, 2006):
488                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
489                            ( 1.0_wp - tau_cloud**0.68_wp )**3
490!
491!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
492!--                (Use constant nu_c = 1.0_wp instead?)
493                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
494!
495!--                Mean weight of cloud droplets:
496                   xc = hyrho(k) * qc(k,j,i) / nc_const
497!
498!--                Parameterized turbulence effects on autoconversion (Seifert,
499!--                Nuijens and Stevens, 2010)
500                   IF ( collision_turbulence )  THEN
501!
502!--                   Weight averaged radius of cloud droplets:
503                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
504
505                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
506                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
507                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
508!
509!--                   Mixing length (neglecting distance to ground and
510!--                   stratification)
511                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
512!
513!--                   Limit dissipation rate according to Seifert, Nuijens and
514!--                   Stevens (2010)
515                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
516!
517!--                   Compute Taylor-microscale Reynolds number:
518                      re_lambda = 6.0_wp / 11.0_wp *                           &
519                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
520                                  SQRT( 15.0_wp / kin_vis_air ) *              &
521                                  dissipation**( 1.0_wp / 6.0_wp )
522!
523!--                   The factor of 1.0E4 is needed to convert the dissipation
524!--                   rate from m2 s-3 to cm2 s-3.
525                      k_au = k_au * ( 1.0_wp +                                 &
526                                      dissipation * 1.0E4_wp *                 &
527                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
528                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &     
529                                                                      r_cc ) / &
530                                                        sigma_cc )**2          &
531                                                      ) + beta_cc              &
532                                      )                                        &
533                                    )
534                   ENDIF
535!
536!--                Autoconversion rate (Seifert and Beheng, 2006):
537                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
538                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
539                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
540                             rho_surface
541                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
542
543                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
544                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
545                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
546
547                ENDIF
548
549             ENDDO
550          ENDDO
551       ENDDO
552
553       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
554
555    END SUBROUTINE autoconversion
556
557
558!------------------------------------------------------------------------------!
559! Description:
560! ------------
561!> Autoconversion process (Kessler, 1969).
562!------------------------------------------------------------------------------!
563    SUBROUTINE autoconversion_kessler
564
565       USE arrays_3d,                                                          &
566           ONLY:  dzw, pt, prr, q, qc
567
568       USE cloud_parameters,                                                   &
569           ONLY:  l_d_cp, pt_d_t
570
571       USE indices,                                                            &
572           ONLY:  nxl, nxr, nyn, nys, nzb_s_inner, nzt
573
574       USE kinds
575
576
577       IMPLICIT NONE
578
579       INTEGER(iwp) ::  i !<
580       INTEGER(iwp) ::  j !<
581       INTEGER(iwp) ::  k !<
582
583       REAL(wp)    ::  dqdt_precip !<
584
585       DO  i = nxl, nxr
586          DO  j = nys, nyn
587             DO  k = nzb_s_inner(j,i)+1, nzt
588
589                IF ( qc(k,j,i) > ql_crit )  THEN
590                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
591                ELSE
592                   dqdt_precip = 0.0_wp
593                ENDIF
594
595                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
596                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro
597                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
598                                        pt_d_t(k)
599
600!
601!--             Compute the rain rate (stored on surface grid point)
602                prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +        &
603                                            dqdt_precip * dzw(k)
604
605             ENDDO
606          ENDDO
607       ENDDO
608
609   END SUBROUTINE autoconversion_kessler
610
611
612!------------------------------------------------------------------------------!
613! Description:
614! ------------
615!> Accretion rate (Seifert and Beheng, 2006).
616!------------------------------------------------------------------------------!
617    SUBROUTINE accretion
618
619       USE arrays_3d,                                                          &
620           ONLY:  diss, qc, qr
621
622       USE cloud_parameters,                                                   &
623           ONLY:  hyrho
624
625       USE control_parameters,                                                 &
626           ONLY:  rho_surface
627
628       USE cpulog,                                                             &
629           ONLY:  cpu_log, log_point_s
630
631       USE indices,                                                            &
632           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
633
634       USE kinds
635
636       IMPLICIT NONE
637
638       INTEGER(iwp) ::  i                 !<
639       INTEGER(iwp) ::  j                 !<
640       INTEGER(iwp) ::  k                 !<
641
642       REAL(wp)     ::  accr              !<
643       REAL(wp)     ::  k_cr              !<
644       REAL(wp)     ::  phi_ac            !<
645       REAL(wp)     ::  tau_cloud         !<
646
647       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
648
649       DO  i = nxl, nxr
650          DO  j = nys, nyn
651             DO  k = nzb_s_inner(j,i)+1, nzt
652
653                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
654!
655!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
656                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 
657!
658!--                Universal function for accretion process (Seifert and
659!--                Beheng, 2001):
660                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
661!
662!--                Parameterized turbulence effects on autoconversion (Seifert,
663!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
664!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
665                   IF ( collision_turbulence )  THEN
666                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
667                                       MIN( 600.0_wp,                          &
668                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
669                                     )
670                   ELSE
671                      k_cr = k_cr0                       
672                   ENDIF
673!
674!--                Accretion rate (Seifert and Beheng, 2006):
675                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
676                          SQRT( rho_surface * hyrho(k) )
677                   accr = MIN( accr, qc(k,j,i) / dt_micro )
678
679                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro
680                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
681
682                ENDIF
683
684             ENDDO
685          ENDDO
686       ENDDO
687
688       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
689
690    END SUBROUTINE accretion
691
692
693!------------------------------------------------------------------------------!
694! Description:
695! ------------
696!> Collisional breakup rate (Seifert, 2008).
697!------------------------------------------------------------------------------!
698    SUBROUTINE selfcollection_breakup
699
700       USE arrays_3d,                                                          &
701           ONLY:  nr, qr
702
703       USE cloud_parameters,                                                   &
704           ONLY:  hyrho
705
706       USE control_parameters,                                                 &
707           ONLY:  rho_surface
708
709       USE cpulog,                                                             &
710           ONLY:  cpu_log, log_point_s
711
712       USE indices,                                                            &
713           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
714
715       USE kinds
716   
717       IMPLICIT NONE
718
719       INTEGER(iwp) ::  i                 !<
720       INTEGER(iwp) ::  j                 !<
721       INTEGER(iwp) ::  k                 !<
722
723       REAL(wp)     ::  breakup           !<
724       REAL(wp)     ::  dr                !<
725       REAL(wp)     ::  phi_br            !<
726       REAL(wp)     ::  selfcoll          !<
727
728       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
729
730       DO  i = nxl, nxr
731          DO  j = nys, nyn
732             DO  k = nzb_s_inner(j,i)+1, nzt
733                IF ( qr(k,j,i) > eps_sb )  THEN
734!
735!--                Selfcollection rate (Seifert and Beheng, 2001):
736                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
737                              SQRT( hyrho(k) * rho_surface )
738!
739!--                Weight averaged diameter of rain drops:
740                   dr = ( hyrho(k) * qr(k,j,i) /                               &
741                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
742!
743!--                Collisional breakup rate (Seifert, 2008):
744                   IF ( dr >= 0.3E-3_wp )  THEN
745                      phi_br  = k_br * ( dr - 1.1E-3_wp )
746                      breakup = selfcoll * ( phi_br + 1.0_wp )
747                   ELSE
748                      breakup = 0.0_wp
749                   ENDIF
750
751                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
752                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
753
754                ENDIF         
755             ENDDO
756          ENDDO
757       ENDDO
758
759       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
760
761    END SUBROUTINE selfcollection_breakup
762
763
764!------------------------------------------------------------------------------!
765! Description:
766! ------------
767!> Evaporation of precipitable water. Condensation is neglected for
768!> precipitable water.
769!------------------------------------------------------------------------------!
770    SUBROUTINE evaporation_rain
771
772       USE arrays_3d,                                                          &
773           ONLY:  hyp, nr, pt, q,  qc, qr
774
775       USE cloud_parameters,                                                   &
776           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
777
778       USE constants,                                                          &
779           ONLY:  pi
780
781       USE cpulog,                                                             &
782           ONLY:  cpu_log, log_point_s
783
784       USE indices,                                                            &
785           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
786
787       USE kinds
788
789       IMPLICIT NONE
790
791       INTEGER(iwp) ::  i                 !<
792       INTEGER(iwp) ::  j                 !<
793       INTEGER(iwp) ::  k                 !<
794
795       REAL(wp)     ::  alpha             !<
796       REAL(wp)     ::  dr                !<
797       REAL(wp)     ::  e_s               !<
798       REAL(wp)     ::  evap              !<
799       REAL(wp)     ::  evap_nr           !<
800       REAL(wp)     ::  f_vent            !<
801       REAL(wp)     ::  g_evap            !<
802       REAL(wp)     ::  lambda_r          !<
803       REAL(wp)     ::  mu_r              !<
804       REAL(wp)     ::  mu_r_2            !<
805       REAL(wp)     ::  mu_r_5d2          !<
806       REAL(wp)     ::  nr_0              !<
807       REAL(wp)     ::  q_s               !<
808       REAL(wp)     ::  sat               !<
809       REAL(wp)     ::  t_l               !<
810       REAL(wp)     ::  temp              !<
811       REAL(wp)     ::  xr                !<
812
813       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
814
815       DO  i = nxl, nxr
816          DO  j = nys, nyn
817             DO  k = nzb_s_inner(j,i)+1, nzt
818                IF ( qr(k,j,i) > eps_sb )  THEN
819!
820!--                Actual liquid water temperature:
821                   t_l = t_d_pt(k) * pt(k,j,i)
822!
823!--                Saturation vapor pressure at t_l:
824                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
825                                          ( t_l - 35.86_wp )                   &
826                                        )
827!
828!--                Computation of saturation humidity:
829                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
830                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
831                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
832                           ( 1.0_wp + alpha * q_s )
833!
834!--                Supersaturation:
835                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
836!
837!--                Evaporation needs only to be calculated in subsaturated regions
838                   IF ( sat < 0.0_wp )  THEN
839!
840!--                   Actual temperature:
841                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
842             
843                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
844                                          l_v / ( thermal_conductivity_l * temp ) &
845                                          + r_v * temp / ( diff_coeff_l * e_s )   &
846                                        )
847!
848!--                   Mean weight of rain drops
849                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
850!
851!--                   Weight averaged diameter of rain drops:
852                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
853!
854!--                   Compute ventilation factor and intercept parameter
855!--                   (Seifert and Beheng, 2006; Seifert, 2008):
856                      IF ( ventilation_effect )  THEN
857!
858!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
859!--                      2005; Stevens and Seifert, 2008):
860                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
861                                                          ( dr - 1.4E-3_wp ) ) )
862!
863!--                      Slope parameter of gamma distribution (Seifert, 2008):
864                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
865                                      ( mu_r + 1.0_wp )                        &
866                                    )**( 1.0_wp / 3.0_wp ) / dr
867
868                         mu_r_2   = mu_r + 2.0_wp
869                         mu_r_5d2 = mu_r + 2.5_wp 
870
871                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
872                                  lambda_r**( -mu_r_2 ) + b_vent *              &
873                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
874                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
875                                  ( 1.0_wp -                                    &
876                                    0.5_wp * ( b_term / a_term ) *              &
877                                    ( lambda_r / ( c_term + lambda_r )          &
878                                    )**mu_r_5d2 -                               &
879                                    0.125_wp * ( b_term / a_term )**2 *         &
880                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
881                                    )**mu_r_5d2 -                               &
882                                    0.0625_wp * ( b_term / a_term )**3 *        &
883                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
884                                    )**mu_r_5d2 -                               &
885                                    0.0390625_wp * ( b_term / a_term )**4 *     & 
886                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
887                                    )**mu_r_5d2                                 &
888                                  )
889
890                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
891                                  gamm( mu_r + 1.0_wp ) 
892                      ELSE
893                         f_vent = 1.0_wp
894                         nr_0   = nr(k,j,i) * dr
895                      ENDIF
896   !
897   !--                Evaporation rate of rain water content (Seifert and
898   !--                Beheng, 2006):
899                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
900                                hyrho(k)
901                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
902                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
903                                     -nr(k,j,i) / dt_micro )
904
905                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
906                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
907
908                   ENDIF
909                ENDIF         
910
911             ENDDO
912          ENDDO
913       ENDDO
914
915       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
916
917    END SUBROUTINE evaporation_rain
918
919
920!------------------------------------------------------------------------------!
921! Description:
922! ------------
923!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
924!------------------------------------------------------------------------------!
925    SUBROUTINE sedimentation_cloud
926
927       USE arrays_3d,                                                          &
928           ONLY:  ddzu, dzu, pt, prr, q, qc
929
930       USE cloud_parameters,                                                   &
931           ONLY:  hyrho, l_d_cp, pt_d_t
932
933       USE control_parameters,                                                 &
934           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
935
936       USE cpulog,                                                             &
937           ONLY:  cpu_log, log_point_s
938
939       USE indices,                                                            &
940           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
941
942       USE kinds
943
944       USE statistics,                                                         &
945           ONLY:  weight_substep
946   
947
948       IMPLICIT NONE
949
950       INTEGER(iwp) ::  i             !<
951       INTEGER(iwp) ::  j             !<
952       INTEGER(iwp) ::  k             !<
953
954       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
955
956       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
957
958       sed_qc(nzt+1) = 0.0_wp
959
960       DO  i = nxl, nxr
961          DO  j = nys, nyn
962             DO  k = nzt, nzb_s_inner(j,i)+1, -1
963
964                IF ( qc(k,j,i) > eps_sb )  THEN
965                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
966                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
967                ELSE
968                   sed_qc(k) = 0.0_wp
969                ENDIF
970
971                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
972                                            dt_micro + sed_qc(k+1)             &
973                               )
974
975                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
976                                        ddzu(k+1) / hyrho(k) * dt_micro
977                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
978                                        ddzu(k+1) / hyrho(k) * dt_micro
979                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
980                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
981                                        pt_d_t(k) * dt_micro
982
983!
984!--             Compute the precipitation rate due to cloud (fog) droplets
985                IF ( call_microphysics_at_all_substeps )  THEN
986                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
987                                * weight_substep(intermediate_timestep_count)
988                ELSE
989                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
990                ENDIF
991
992             ENDDO
993          ENDDO
994       ENDDO
995
996       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
997
998    END SUBROUTINE sedimentation_cloud
999
1000
1001!------------------------------------------------------------------------------!
1002! Description:
1003! ------------
1004!> Computation of sedimentation flux. Implementation according to Stevens
1005!> and Seifert (2008). Code is based on UCLA-LES.
1006!------------------------------------------------------------------------------!
1007    SUBROUTINE sedimentation_rain
1008
1009       USE arrays_3d,                                                          &
1010           ONLY:  ddzu, dzu, nr, pt, prr, q, qr
1011
1012       USE cloud_parameters,                                                   &
1013           ONLY:  hyrho, l_d_cp, pt_d_t
1014
1015       USE control_parameters,                                                 &
1016           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
1017       USE cpulog,                                                             &
1018           ONLY:  cpu_log, log_point_s
1019
1020       USE indices,                                                            &
1021           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
1022
1023       USE kinds
1024
1025       USE statistics,                                                         &
1026           ONLY:  weight_substep
1027       
1028       IMPLICIT NONE
1029
1030       INTEGER(iwp) ::  i                          !<
1031       INTEGER(iwp) ::  j                          !<
1032       INTEGER(iwp) ::  k                          !<
1033       INTEGER(iwp) ::  k_run                      !<
1034
1035       REAL(wp)     ::  c_run                      !<
1036       REAL(wp)     ::  d_max                      !<
1037       REAL(wp)     ::  d_mean                     !<
1038       REAL(wp)     ::  d_min                      !<
1039       REAL(wp)     ::  dr                         !<
1040       REAL(wp)     ::  flux                       !<
1041       REAL(wp)     ::  lambda_r                   !<
1042       REAL(wp)     ::  mu_r                       !<
1043       REAL(wp)     ::  z_run                      !<
1044
1045       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1046       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1047       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1048       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1049       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1050       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1051       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1052       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
1053
1054       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
1055
1056!
1057!--    Compute velocities
1058       DO  i = nxl, nxr
1059          DO  j = nys, nyn
1060             DO  k = nzb_s_inner(j,i)+1, nzt
1061                IF ( qr(k,j,i) > eps_sb )  THEN
1062!
1063!--                Weight averaged diameter of rain drops:
1064                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1065                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1066!
1067!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1068!--                Stevens and Seifert, 2008):
1069                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
1070                                                     ( dr - 1.4E-3_wp ) ) )
1071!
1072!--                Slope parameter of gamma distribution (Seifert, 2008):
1073                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1074                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1075
1076                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1077                                               a_term - b_term * ( 1.0_wp +    &
1078                                                  c_term /                     &
1079                                                  lambda_r )**( -1.0_wp *      &
1080                                                  ( mu_r + 1.0_wp ) )          &
1081                                              )                                &
1082                                )
1083
1084                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1085                                               a_term - b_term * ( 1.0_wp +    &
1086                                                  c_term /                     &
1087                                                  lambda_r )**( -1.0_wp *      &
1088                                                  ( mu_r + 4.0_wp ) )          &
1089                                             )                                 &
1090                                )
1091                ELSE
1092                   w_nr(k) = 0.0_wp
1093                   w_qr(k) = 0.0_wp
1094                ENDIF
1095             ENDDO
1096!
1097!--          Adjust boundary values
1098             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
1099             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
1100             w_nr(nzt+1) = 0.0_wp
1101             w_qr(nzt+1) = 0.0_wp
1102!
1103!--          Compute Courant number
1104             DO  k = nzb_s_inner(j,i)+1, nzt
1105                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
1106                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
1107                          dt_micro * ddzu(k)
1108                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
1109                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
1110                          dt_micro * ddzu(k)
1111             ENDDO     
1112!
1113!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1114             IF ( limiter_sedimentation )  THEN
1115
1116                DO k = nzb_s_inner(j,i)+1, nzt
1117                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
1118                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1119                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1120
1121                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1122                                                              2.0_wp * d_max,  &
1123                                                              ABS( d_mean ) )
1124
1125                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
1126                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1127                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1128
1129                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1130                                                              2.0_wp * d_max,  &
1131                                                              ABS( d_mean ) )
1132                ENDDO
1133
1134             ELSE
1135
1136                nr_slope = 0.0_wp
1137                qr_slope = 0.0_wp
1138
1139             ENDIF
1140
1141             sed_nr(nzt+1) = 0.0_wp
1142             sed_qr(nzt+1) = 0.0_wp
1143!
1144!--          Compute sedimentation flux
1145             DO  k = nzt, nzb_s_inner(j,i)+1, -1
1146!
1147!--             Sum up all rain drop number densities which contribute to the flux
1148!--             through k-1/2
1149                flux  = 0.0_wp
1150                z_run = 0.0_wp ! height above z(k)
1151                k_run = k
1152                c_run = MIN( 1.0_wp, c_nr(k) )
1153                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1154                   flux  = flux + hyrho(k_run) *                               &
1155                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
1156                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
1157                   z_run = z_run + dzu(k_run)
1158                   k_run = k_run + 1
1159                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
1160                ENDDO
1161!
1162!--             It is not allowed to sediment more rain drop number density than
1163!--             available
1164                flux = MIN( flux,                                              &
1165                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
1166                            dt_micro                                           &
1167                          )
1168
1169                sed_nr(k) = flux / dt_micro
1170                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
1171                                        ddzu(k+1) / hyrho(k) * dt_micro
1172!
1173!--             Sum up all rain water content which contributes to the flux
1174!--             through k-1/2
1175                flux  = 0.0_wp
1176                z_run = 0.0_wp ! height above z(k)
1177                k_run = k
1178                c_run = MIN( 1.0_wp, c_qr(k) )
1179
1180                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1181
1182                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
1183                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
1184                                  0.5_wp ) * c_run * dzu(k_run)
1185                   z_run = z_run + dzu(k_run)
1186                   k_run = k_run + 1
1187                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
1188
1189                ENDDO
1190!
1191!--             It is not allowed to sediment more rain water content than
1192!--             available
1193                flux = MIN( flux,                                              &
1194                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
1195                            dt_micro                                           &
1196                          )
1197
1198                sed_qr(k) = flux / dt_micro
1199
1200                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
1201                                        ddzu(k+1) / hyrho(k) * dt_micro
1202                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
1203                                        ddzu(k+1) / hyrho(k) * dt_micro
1204                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
1205                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1206                                        pt_d_t(k) * dt_micro
1207!
1208!--             Compute the rain rate
1209                IF ( call_microphysics_at_all_substeps )  THEN
1210                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
1211                                * weight_substep(intermediate_timestep_count)
1212                ELSE
1213                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
1214                ENDIF
1215
1216             ENDDO
1217          ENDDO
1218       ENDDO
1219
1220       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1221
1222    END SUBROUTINE sedimentation_rain
1223
1224
1225!------------------------------------------------------------------------------!
1226! Description:
1227! ------------
1228!> Computation of the precipitation amount due to gravitational settling of
1229!> rain and cloud (fog) droplets
1230!------------------------------------------------------------------------------!
1231    SUBROUTINE calc_precipitation_amount
1232
1233       USE arrays_3d,                                                          &
1234           ONLY:  precipitation_amount, prr
1235
1236       USE cloud_parameters,                                                   &
1237           ONLY:  hyrho
1238
1239       USE control_parameters,                                                 &
1240           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
1241                  intermediate_timestep_count, intermediate_timestep_count_max,&
1242                  precipitation_amount_interval, time_do2d_xy
1243
1244       USE indices,                                                            &
1245           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner
1246
1247       USE kinds
1248
1249       IMPLICIT NONE
1250
1251       INTEGER(iwp) ::  i                          !:
1252       INTEGER(iwp) ::  j                          !:
1253
1254
1255       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1256            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
1257            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1258       THEN
1259
1260          DO  i = nxl, nxr
1261             DO  j = nys, nyn
1262
1263                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
1264                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
1265                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
1266
1267             ENDDO
1268          ENDDO
1269       ENDIF
1270
1271    END SUBROUTINE calc_precipitation_amount
1272
1273
1274!------------------------------------------------------------------------------!
1275! Description:
1276! ------------
1277!> Control of microphysics for grid points i,j
1278!------------------------------------------------------------------------------!
1279
1280    SUBROUTINE microphysics_control_ij( i, j )
1281
1282       USE arrays_3d,                                                          &
1283           ONLY:  hyp, nr, pt, pt_init, prr, q, qc, qr, zu
1284
1285       USE cloud_parameters,                                                   &
1286           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
1287
1288       USE control_parameters,                                                 &
1289           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
1290                  intermediate_timestep_count, large_scale_forcing,            &
1291                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
1292                  pt_surface, rho_surface, surface_pressure
1293
1294       USE indices,                                                            &
1295           ONLY:  nzb, nzt
1296
1297       USE kinds
1298
1299       USE statistics,                                                         &
1300           ONLY:  weight_pres
1301
1302       IMPLICIT NONE
1303
1304       INTEGER(iwp) ::  i                 !<
1305       INTEGER(iwp) ::  j                 !<
1306       INTEGER(iwp) ::  k                 !<
1307
1308       REAL(wp)     ::  t_surface         !<
1309
1310       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
1311!
1312!--       Calculate:
1313!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1314!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1315!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1316          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1317          DO  k = nzb, nzt+1
1318             hyp(k)    = surface_pressure * 100.0_wp * &
1319                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1320             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1321             t_d_pt(k) = 1.0_wp / pt_d_t(k)
1322             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
1323          ENDDO
1324!
1325!--       Compute reference density
1326          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1327       ENDIF
1328
1329!
1330!--    Compute length of time step
1331       IF ( call_microphysics_at_all_substeps )  THEN
1332          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1333       ELSE
1334          dt_micro = dt_3d
1335       ENDIF
1336
1337!
1338!--    Use 1d arrays
1339       q_1d(:)  = q(:,j,i)
1340       pt_1d(:) = pt(:,j,i)
1341       qc_1d(:) = qc(:,j,i)
1342       nc_1d(:) = nc_const
1343       IF ( microphysics_seifert )  THEN
1344          qr_1d(:) = qr(:,j,i)
1345          nr_1d(:) = nr(:,j,i)
1346       ENDIF
1347
1348!
1349!--    Reset precipitation rate
1350       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1351
1352!
1353!--    Compute cloud physics
1354       IF( microphysics_kessler )  THEN
1355
1356          CALL autoconversion_kessler( i,j )
1357          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1358
1359       ELSEIF ( microphysics_seifert )  THEN
1360
1361          CALL adjust_cloud( i,j )
1362          CALL autoconversion( i,j )
1363          CALL accretion( i,j )
1364          CALL selfcollection_breakup( i,j )
1365          CALL evaporation_rain( i,j )
1366          CALL sedimentation_rain( i,j )
1367          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1368
1369       ENDIF
1370
1371       CALL calc_precipitation_amount( i,j )
1372
1373!
1374!--    Store results on the 3d arrays
1375       q(:,j,i)  = q_1d(:)
1376       pt(:,j,i) = pt_1d(:)
1377       IF ( microphysics_seifert )  THEN
1378          qr(:,j,i) = qr_1d(:)
1379          nr(:,j,i) = nr_1d(:)
1380       ENDIF
1381
1382    END SUBROUTINE microphysics_control_ij
1383
1384!------------------------------------------------------------------------------!
1385! Description:
1386! ------------
1387!> Adjust number of raindrops to avoid nonlinear effects in
1388!> sedimentation and evaporation of rain drops due to too small or
1389!> too big weights of rain drops (Stevens and Seifert, 2008).
1390!> The same procedure is applied to cloud droplets if they are determined
1391!> prognostically. Call for grid point i,j
1392!------------------------------------------------------------------------------!
1393    SUBROUTINE adjust_cloud_ij( i, j )
1394
1395       USE cloud_parameters,                                                   &
1396           ONLY:  hyrho
1397
1398       USE indices,                                                            &
1399           ONLY:  nzb_s_inner, nzt
1400
1401       USE kinds
1402
1403       IMPLICIT NONE
1404
1405       INTEGER(iwp) ::  i                 !<
1406       INTEGER(iwp) ::  j                 !<
1407       INTEGER(iwp) ::  k                 !<
1408
1409       DO  k = nzb_s_inner(j,i)+1, nzt
1410
1411          IF ( qr_1d(k) <= eps_sb )  THEN
1412             qr_1d(k) = 0.0_wp
1413             nr_1d(k) = 0.0_wp
1414          ELSE
1415!
1416!--          Adjust number of raindrops to avoid nonlinear effects in
1417!--          sedimentation and evaporation of rain drops due to too small or
1418!--          too big weights of rain drops (Stevens and Seifert, 2008).
1419             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
1420                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
1421             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
1422                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
1423             ENDIF
1424
1425          ENDIF
1426
1427       ENDDO
1428
1429    END SUBROUTINE adjust_cloud_ij
1430
1431
1432!------------------------------------------------------------------------------!
1433! Description:
1434! ------------
1435!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
1436!------------------------------------------------------------------------------!
1437    SUBROUTINE autoconversion_ij( i, j )
1438
1439       USE arrays_3d,                                                          &
1440           ONLY:  diss, dzu
1441
1442       USE cloud_parameters,                                                   &
1443           ONLY:  hyrho
1444
1445       USE control_parameters,                                                 &
1446           ONLY:  rho_surface
1447
1448       USE grid_variables,                                                     &
1449           ONLY:  dx, dy
1450
1451       USE indices,                                                            &
1452           ONLY:  nzb_s_inner, nzt
1453
1454       USE kinds
1455
1456       IMPLICIT NONE
1457
1458       INTEGER(iwp) ::  i                 !<
1459       INTEGER(iwp) ::  j                 !<
1460       INTEGER(iwp) ::  k                 !<
1461
1462       REAL(wp)     ::  alpha_cc          !<                   
1463       REAL(wp)     ::  autocon           !<
1464       REAL(wp)     ::  dissipation       !<
1465       REAL(wp)     ::  k_au              !<
1466       REAL(wp)     ::  l_mix             !<
1467       REAL(wp)     ::  nu_c              !<
1468       REAL(wp)     ::  phi_au            !<
1469       REAL(wp)     ::  r_cc              !<
1470       REAL(wp)     ::  rc                !<
1471       REAL(wp)     ::  re_lambda         !<
1472       REAL(wp)     ::  sigma_cc          !<
1473       REAL(wp)     ::  tau_cloud         !<
1474       REAL(wp)     ::  xc                !<
1475
1476       DO  k = nzb_s_inner(j,i)+1, nzt
1477
1478          IF ( qc_1d(k) > eps_sb )  THEN
1479
1480             k_au = k_cc / ( 20.0_wp * x0 )
1481!
1482!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1483!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1484             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
1485!
1486!--          Universal function for autoconversion process
1487!--          (Seifert and Beheng, 2006):
1488             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
1489!
1490!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
1491!--          (Use constant nu_c = 1.0_wp instead?)
1492             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
1493!
1494!--          Mean weight of cloud droplets:
1495             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
1496!
1497!--          Parameterized turbulence effects on autoconversion (Seifert,
1498!--          Nuijens and Stevens, 2010)
1499             IF ( collision_turbulence )  THEN
1500!
1501!--             Weight averaged radius of cloud droplets:
1502                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
1503
1504                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1505                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1506                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
1507!
1508!--             Mixing length (neglecting distance to ground and stratification)
1509                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
1510!
1511!--             Limit dissipation rate according to Seifert, Nuijens and
1512!--             Stevens (2010)
1513                dissipation = MIN( 0.06_wp, diss(k,j,i) )
1514!
1515!--             Compute Taylor-microscale Reynolds number:
1516                re_lambda = 6.0_wp / 11.0_wp *                                 &
1517                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
1518                            SQRT( 15.0_wp / kin_vis_air ) *                    &
1519                            dissipation**( 1.0_wp / 6.0_wp )
1520!
1521!--             The factor of 1.0E4 is needed to convert the dissipation rate
1522!--             from m2 s-3 to cm2 s-3.
1523                k_au = k_au * ( 1.0_wp +                                       &
1524                                dissipation * 1.0E4_wp *                       &
1525                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
1526                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
1527                                                  sigma_cc )**2                &
1528                                                ) + beta_cc                    &
1529                                )                                              &
1530                              )
1531             ENDIF
1532!
1533!--          Autoconversion rate (Seifert and Beheng, 2006):
1534             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
1535                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
1536                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
1537                       rho_surface
1538             autocon = MIN( autocon, qc_1d(k) / dt_micro )
1539
1540             qr_1d(k) = qr_1d(k) + autocon * dt_micro
1541             qc_1d(k) = qc_1d(k) - autocon * dt_micro
1542             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
1543
1544          ENDIF
1545
1546       ENDDO
1547
1548    END SUBROUTINE autoconversion_ij
1549
1550!------------------------------------------------------------------------------!
1551! Description:
1552! ------------
1553!> Autoconversion process (Kessler, 1969).
1554!------------------------------------------------------------------------------!
1555    SUBROUTINE autoconversion_kessler_ij( i, j )
1556
1557       USE arrays_3d,                                                          &
1558           ONLY:  dzw, prr
1559
1560       USE cloud_parameters,                                                   &
1561           ONLY:  l_d_cp, pt_d_t
1562
1563       USE indices,                                                            &
1564           ONLY:  nzb_s_inner, nzt
1565
1566       USE kinds
1567
1568
1569       IMPLICIT NONE
1570
1571       INTEGER(iwp) ::  i !<
1572       INTEGER(iwp) ::  j !<
1573       INTEGER(iwp) ::  k !<
1574
1575       REAL(wp)    ::  dqdt_precip !<
1576
1577       DO  k = nzb_s_inner(j,i)+1, nzt
1578
1579          IF ( qc_1d(k) > ql_crit )  THEN
1580             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
1581          ELSE
1582             dqdt_precip = 0.0_wp
1583          ENDIF
1584
1585          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
1586          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro
1587          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
1588
1589!
1590!--       Compute the rain rate (stored on surface grid point)
1591          prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +              &
1592                                      dqdt_precip * dzw(k)
1593
1594       ENDDO
1595
1596    END SUBROUTINE autoconversion_kessler_ij
1597
1598!------------------------------------------------------------------------------!
1599! Description:
1600! ------------
1601!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
1602!------------------------------------------------------------------------------!
1603    SUBROUTINE accretion_ij( i, j )
1604
1605       USE arrays_3d,                                                          &
1606           ONLY:  diss
1607
1608       USE cloud_parameters,                                                   &
1609           ONLY:  hyrho
1610
1611       USE control_parameters,                                                 &
1612           ONLY:  rho_surface
1613
1614       USE indices,                                                            &
1615           ONLY:  nzb_s_inner, nzt
1616
1617       USE kinds
1618
1619       IMPLICIT NONE
1620
1621       INTEGER(iwp) ::  i                 !<
1622       INTEGER(iwp) ::  j                 !<
1623       INTEGER(iwp) ::  k                 !<
1624
1625       REAL(wp)     ::  accr              !<
1626       REAL(wp)     ::  k_cr              !<
1627       REAL(wp)     ::  phi_ac            !<
1628       REAL(wp)     ::  tau_cloud         !<
1629
1630       DO  k = nzb_s_inner(j,i)+1, nzt
1631          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
1632!
1633!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1634             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
1635!
1636!--          Universal function for accretion process
1637!--          (Seifert and Beheng, 2001):
1638             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1639!
1640!--          Parameterized turbulence effects on autoconversion (Seifert,
1641!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1642!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1643             IF ( collision_turbulence )  THEN
1644                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
1645                                 MIN( 600.0_wp,                                &
1646                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
1647                               )
1648             ELSE
1649                k_cr = k_cr0                       
1650             ENDIF
1651!
1652!--          Accretion rate (Seifert and Beheng, 2006):
1653             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
1654             accr = MIN( accr, qc_1d(k) / dt_micro )
1655
1656             qr_1d(k) = qr_1d(k) + accr * dt_micro
1657             qc_1d(k) = qc_1d(k) - accr * dt_micro
1658
1659          ENDIF
1660
1661       ENDDO
1662
1663    END SUBROUTINE accretion_ij
1664
1665
1666!------------------------------------------------------------------------------!
1667! Description:
1668! ------------
1669!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
1670!------------------------------------------------------------------------------!
1671    SUBROUTINE selfcollection_breakup_ij( i, j )
1672
1673       USE cloud_parameters,                                                   &
1674           ONLY:  hyrho
1675
1676       USE control_parameters,                                                 &
1677           ONLY:  rho_surface
1678
1679       USE indices,                                                            &
1680           ONLY:  nzb_s_inner, nzt
1681
1682       USE kinds
1683   
1684       IMPLICIT NONE
1685
1686       INTEGER(iwp) ::  i                 !<
1687       INTEGER(iwp) ::  j                 !<
1688       INTEGER(iwp) ::  k                 !<
1689
1690       REAL(wp)     ::  breakup           !<
1691       REAL(wp)     ::  dr                !<
1692       REAL(wp)     ::  phi_br            !<
1693       REAL(wp)     ::  selfcoll          !<
1694
1695       DO  k = nzb_s_inner(j,i)+1, nzt
1696          IF ( qr_1d(k) > eps_sb )  THEN
1697!
1698!--          Selfcollection rate (Seifert and Beheng, 2001):
1699             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
1700!
1701!--          Weight averaged diameter of rain drops:
1702             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1703!
1704!--          Collisional breakup rate (Seifert, 2008):
1705             IF ( dr >= 0.3E-3_wp )  THEN
1706                phi_br  = k_br * ( dr - 1.1E-3_wp )
1707                breakup = selfcoll * ( phi_br + 1.0_wp )
1708             ELSE
1709                breakup = 0.0_wp
1710             ENDIF
1711
1712             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1713             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
1714
1715          ENDIF         
1716       ENDDO
1717
1718    END SUBROUTINE selfcollection_breakup_ij
1719
1720
1721!------------------------------------------------------------------------------!
1722! Description:
1723! ------------
1724!> Evaporation of precipitable water. Condensation is neglected for
1725!> precipitable water. Call for grid point i,j
1726!------------------------------------------------------------------------------!
1727    SUBROUTINE evaporation_rain_ij( i, j )
1728
1729       USE arrays_3d,                                                          &
1730           ONLY:  hyp
1731
1732       USE cloud_parameters,                                                   &
1733           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
1734
1735       USE constants,                                                          &
1736           ONLY:  pi
1737
1738       USE indices,                                                            &
1739           ONLY:  nzb_s_inner, nzt
1740
1741       USE kinds
1742
1743       IMPLICIT NONE
1744
1745       INTEGER(iwp) ::  i                 !<
1746       INTEGER(iwp) ::  j                 !<
1747       INTEGER(iwp) ::  k                 !<
1748
1749       REAL(wp)     ::  alpha             !<
1750       REAL(wp)     ::  dr                !<
1751       REAL(wp)     ::  e_s               !<
1752       REAL(wp)     ::  evap              !<
1753       REAL(wp)     ::  evap_nr           !<
1754       REAL(wp)     ::  f_vent            !<
1755       REAL(wp)     ::  g_evap            !<
1756       REAL(wp)     ::  lambda_r          !<
1757       REAL(wp)     ::  mu_r              !<
1758       REAL(wp)     ::  mu_r_2            !<
1759       REAL(wp)     ::  mu_r_5d2          !<
1760       REAL(wp)     ::  nr_0              !<
1761       REAL(wp)     ::  q_s               !<
1762       REAL(wp)     ::  sat               !<
1763       REAL(wp)     ::  t_l               !<
1764       REAL(wp)     ::  temp              !<
1765       REAL(wp)     ::  xr                !<
1766
1767       DO  k = nzb_s_inner(j,i)+1, nzt
1768          IF ( qr_1d(k) > eps_sb )  THEN
1769!
1770!--          Actual liquid water temperature:
1771             t_l = t_d_pt(k) * pt_1d(k)
1772!
1773!--          Saturation vapor pressure at t_l:
1774             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
1775                                    ( t_l - 35.86_wp )                         &
1776                                  )
1777!
1778!--          Computation of saturation humidity:
1779             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1780             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1781             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
1782!
1783!--          Supersaturation:
1784             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
1785!
1786!--          Evaporation needs only to be calculated in subsaturated regions
1787             IF ( sat < 0.0_wp )  THEN
1788!
1789!--             Actual temperature:
1790                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1791       
1792                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
1793                                    ( thermal_conductivity_l * temp ) +        &
1794                                    r_v * temp / ( diff_coeff_l * e_s )        &
1795                                  )
1796!
1797!--             Mean weight of rain drops
1798                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
1799!
1800!--             Weight averaged diameter of rain drops:
1801                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1802!
1803!--             Compute ventilation factor and intercept parameter
1804!--             (Seifert and Beheng, 2006; Seifert, 2008):
1805                IF ( ventilation_effect )  THEN
1806!
1807!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1808!--                Stevens and Seifert, 2008):
1809                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1810!
1811!--                Slope parameter of gamma distribution (Seifert, 2008):
1812                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1813                                ( mu_r + 1.0_wp )                              &
1814                              )**( 1.0_wp / 3.0_wp ) / dr
1815
1816                   mu_r_2   = mu_r + 2.0_wp
1817                   mu_r_5d2 = mu_r + 2.5_wp 
1818
1819                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
1820                            b_vent * schmidt_p_1d3 *                           &
1821                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
1822                            lambda_r**( -mu_r_5d2 ) *                          &
1823                            ( 1.0_wp -                                         &
1824                              0.5_wp * ( b_term / a_term ) *                   &
1825                              ( lambda_r / ( c_term + lambda_r )               &
1826                              )**mu_r_5d2 -                                    &
1827                              0.125_wp * ( b_term / a_term )**2 *              &
1828                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
1829                              )**mu_r_5d2 -                                    &
1830                              0.0625_wp * ( b_term / a_term )**3 *             &
1831                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
1832                              )**mu_r_5d2 -                                    &
1833                              0.0390625_wp * ( b_term / a_term )**4 *          & 
1834                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
1835                              )**mu_r_5d2                                      &
1836                            )
1837
1838                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
1839                            gamm( mu_r + 1.0_wp ) 
1840                ELSE
1841                   f_vent = 1.0_wp
1842                   nr_0   = nr_1d(k) * dr
1843                ENDIF
1844!
1845!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
1846                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1847                evap    = MAX( evap, -qr_1d(k) / dt_micro )
1848                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
1849                               -nr_1d(k) / dt_micro )
1850
1851                qr_1d(k) = qr_1d(k) + evap    * dt_micro
1852                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
1853
1854             ENDIF
1855          ENDIF         
1856
1857       ENDDO
1858
1859    END SUBROUTINE evaporation_rain_ij
1860
1861
1862!------------------------------------------------------------------------------!
1863! Description:
1864! ------------
1865!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1866!> Call for grid point i,j
1867!------------------------------------------------------------------------------!
1868    SUBROUTINE sedimentation_cloud_ij( i, j )
1869
1870       USE arrays_3d,                                                          &
1871           ONLY:  ddzu, dzu, prr
1872
1873       USE cloud_parameters,                                                   &
1874           ONLY:  hyrho, l_d_cp, pt_d_t
1875
1876       USE control_parameters,                                                 &
1877           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
1878
1879       USE indices,                                                            &
1880           ONLY:  nzb, nzb_s_inner, nzt
1881
1882       USE kinds
1883       
1884       USE statistics,                                                         &
1885           ONLY:  weight_substep
1886
1887       IMPLICIT NONE
1888
1889       INTEGER(iwp) ::  i             !<
1890       INTEGER(iwp) ::  j             !<
1891       INTEGER(iwp) ::  k             !<
1892
1893       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !<
1894
1895       sed_qc(nzt+1) = 0.0_wp
1896
1897       DO  k = nzt, nzb_s_inner(j,i)+1, -1
1898          IF ( qc_1d(k) > eps_sb )  THEN
1899             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
1900                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
1901          ELSE
1902             sed_qc(k) = 0.0_wp
1903          ENDIF
1904
1905          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
1906                                      dt_micro + sed_qc(k+1)                   &
1907                         )
1908
1909          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
1910                                hyrho(k) * dt_micro
1911          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
1912                                hyrho(k) * dt_micro
1913          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
1914                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1915
1916!
1917!--       Compute the precipitation rate of cloud (fog) droplets
1918          IF ( call_microphysics_at_all_substeps )  THEN
1919             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
1920                             weight_substep(intermediate_timestep_count)
1921          ELSE
1922             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
1923          ENDIF
1924
1925       ENDDO
1926
1927    END SUBROUTINE sedimentation_cloud_ij
1928
1929
1930!------------------------------------------------------------------------------!
1931! Description:
1932! ------------
1933!> Computation of sedimentation flux. Implementation according to Stevens
1934!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
1935!------------------------------------------------------------------------------!
1936    SUBROUTINE sedimentation_rain_ij( i, j )
1937
1938       USE arrays_3d,                                                          &
1939           ONLY:  ddzu, dzu, prr
1940
1941       USE cloud_parameters,                                                   &
1942           ONLY:  hyrho, l_d_cp, pt_d_t
1943
1944       USE control_parameters,                                                 &
1945           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
1946
1947       USE indices,                                                            &
1948           ONLY:  nzb, nzb_s_inner, nzt
1949
1950       USE kinds
1951
1952       USE statistics,                                                         &
1953           ONLY:  weight_substep
1954       
1955       IMPLICIT NONE
1956
1957       INTEGER(iwp) ::  i                          !<
1958       INTEGER(iwp) ::  j                          !<
1959       INTEGER(iwp) ::  k                          !<
1960       INTEGER(iwp) ::  k_run                      !<
1961
1962       REAL(wp)     ::  c_run                      !<
1963       REAL(wp)     ::  d_max                      !<
1964       REAL(wp)     ::  d_mean                     !<
1965       REAL(wp)     ::  d_min                      !<
1966       REAL(wp)     ::  dr                         !<
1967       REAL(wp)     ::  flux                       !<
1968       REAL(wp)     ::  lambda_r                   !<
1969       REAL(wp)     ::  mu_r                       !<
1970       REAL(wp)     ::  z_run                      !<
1971
1972       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1973       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1974       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1975       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1976       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1977       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1978       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1979       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
1980
1981!
1982!--    Compute velocities
1983       DO  k = nzb_s_inner(j,i)+1, nzt
1984          IF ( qr_1d(k) > eps_sb )  THEN
1985!
1986!--          Weight averaged diameter of rain drops:
1987             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1988!
1989!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1990!--          Stevens and Seifert, 2008):
1991             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1992!
1993!--          Slope parameter of gamma distribution (Seifert, 2008):
1994             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
1995                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1996
1997             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
1998                                         a_term - b_term * ( 1.0_wp +          &
1999                                            c_term / lambda_r )**( -1.0_wp *   &
2000                                            ( mu_r + 1.0_wp ) )                &
2001                                        )                                      &
2002                          )
2003             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2004                                         a_term - b_term * ( 1.0_wp +          &
2005                                            c_term / lambda_r )**( -1.0_wp *   &
2006                                            ( mu_r + 4.0_wp ) )                &
2007                                       )                                       &
2008                          )
2009          ELSE
2010             w_nr(k) = 0.0_wp
2011             w_qr(k) = 0.0_wp
2012          ENDIF
2013       ENDDO
2014!
2015!--    Adjust boundary values
2016       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
2017       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
2018       w_nr(nzt+1) = 0.0_wp
2019       w_qr(nzt+1) = 0.0_wp
2020!
2021!--    Compute Courant number
2022       DO  k = nzb_s_inner(j,i)+1, nzt
2023          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
2024                    dt_micro * ddzu(k)
2025          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
2026                    dt_micro * ddzu(k)
2027       ENDDO     
2028!
2029!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2030       IF ( limiter_sedimentation )  THEN
2031
2032          DO k = nzb_s_inner(j,i)+1, nzt
2033             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
2034             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2035             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
2036
2037             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2038                                                        2.0_wp * d_max,        &
2039                                                        ABS( d_mean ) )
2040
2041             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
2042             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2043             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
2044
2045             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2046                                                        2.0_wp * d_max,        &
2047                                                        ABS( d_mean ) )
2048          ENDDO
2049
2050       ELSE
2051
2052          nr_slope = 0.0_wp
2053          qr_slope = 0.0_wp
2054
2055       ENDIF
2056
2057       sed_nr(nzt+1) = 0.0_wp
2058       sed_qr(nzt+1) = 0.0_wp
2059!
2060!--    Compute sedimentation flux
2061       DO  k = nzt, nzb_s_inner(j,i)+1, -1
2062!
2063!--       Sum up all rain drop number densities which contribute to the flux
2064!--       through k-1/2
2065          flux  = 0.0_wp
2066          z_run = 0.0_wp ! height above z(k)
2067          k_run = k
2068          c_run = MIN( 1.0_wp, c_nr(k) )
2069          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
2070             flux  = flux + hyrho(k_run) *                                     &
2071                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
2072                     0.5_wp ) * c_run * dzu(k_run)
2073             z_run = z_run + dzu(k_run)
2074             k_run = k_run + 1
2075             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
2076          ENDDO
2077!
2078!--       It is not allowed to sediment more rain drop number density than
2079!--       available
2080          flux = MIN( flux,                                                    &
2081                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
2082
2083          sed_nr(k) = flux / dt_micro
2084          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
2085                                    hyrho(k) * dt_micro
2086!
2087!--       Sum up all rain water content which contributes to the flux
2088!--       through k-1/2
2089          flux  = 0.0_wp
2090          z_run = 0.0_wp ! height above z(k)
2091          k_run = k
2092          c_run = MIN( 1.0_wp, c_qr(k) )
2093
2094          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
2095
2096             flux  = flux + hyrho(k_run) *                                     &
2097                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
2098                     0.5_wp ) * c_run * dzu(k_run)
2099             z_run = z_run + dzu(k_run)
2100             k_run = k_run + 1
2101             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
2102
2103          ENDDO
2104!
2105!--       It is not allowed to sediment more rain water content than available
2106          flux = MIN( flux,                                                    &
2107                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
2108
2109          sed_qr(k) = flux / dt_micro
2110
2111          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2112                                hyrho(k) * dt_micro
2113          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2114                                hyrho(k) * dt_micro
2115          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2116                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
2117!
2118!--       Compute the rain rate
2119          IF ( call_microphysics_at_all_substeps )  THEN
2120             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
2121                          * weight_substep(intermediate_timestep_count)
2122          ELSE
2123             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
2124          ENDIF
2125
2126       ENDDO
2127
2128    END SUBROUTINE sedimentation_rain_ij
2129
2130
2131!------------------------------------------------------------------------------!
2132! Description:
2133! ------------
2134!> This subroutine computes the precipitation amount due to gravitational
2135!> settling of rain and cloud (fog) droplets
2136!------------------------------------------------------------------------------!
2137    SUBROUTINE calc_precipitation_amount_ij( i, j )
2138
2139       USE arrays_3d,                                                          &
2140           ONLY:  precipitation_amount, prr
2141
2142       USE cloud_parameters,                                                   &
2143           ONLY:  hyrho
2144
2145       USE control_parameters,                                                 &
2146           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
2147                  intermediate_timestep_count, intermediate_timestep_count_max,&
2148                  precipitation_amount_interval, time_do2d_xy
2149
2150       USE indices,                                                            &
2151           ONLY:  nzb_s_inner
2152
2153       USE kinds
2154
2155       IMPLICIT NONE
2156
2157       INTEGER(iwp) ::  i                          !:
2158       INTEGER(iwp) ::  j                          !:
2159
2160
2161       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
2162            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
2163            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
2164       THEN
2165
2166          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
2167                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
2168                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
2169       ENDIF
2170
2171    END SUBROUTINE calc_precipitation_amount_ij
2172
2173!------------------------------------------------------------------------------!
2174! Description:
2175! ------------
2176!> This function computes the gamma function (Press et al., 1992).
2177!> The gamma function is needed for the calculation of the evaporation
2178!> of rain drops.
2179!------------------------------------------------------------------------------!
2180    FUNCTION gamm( xx ) 
2181
2182       USE kinds
2183
2184       IMPLICIT NONE
2185
2186       INTEGER(iwp) ::  j            !<
2187
2188       REAL(wp)     ::  gamm         !<
2189       REAL(wp)     ::  ser          !<
2190       REAL(wp)     ::  tmp          !<
2191       REAL(wp)     ::  x_gamm       !<
2192       REAL(wp)     ::  xx           !<
2193       REAL(wp)     ::  y_gamm       !<
2194
2195
2196       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
2197       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
2198                                           -86.50532032941677_wp,      &
2199                                            24.01409824083091_wp,      &
2200                                            -1.231739572450155_wp,     &
2201                                             0.1208650973866179E-2_wp, &
2202                                            -0.5395239384953E-5_wp /)     !<
2203
2204       x_gamm = xx
2205       y_gamm = x_gamm
2206       tmp = x_gamm + 5.5_wp
2207       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
2208       ser = 1.000000000190015_wp
2209
2210       DO  j = 1, 6 
2211          y_gamm = y_gamm + 1.0_wp 
2212          ser    = ser + cof( j ) / y_gamm
2213       ENDDO
2214
2215!
2216!--    Until this point the algorithm computes the logarithm of the gamma
2217!--    function. Hence, the exponential function is used. 
2218!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
2219       gamm = EXP( tmp ) * stp * ser / x_gamm
2220
2221       RETURN
2222 
2223    END FUNCTION gamm
2224
2225 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.