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

Last change on this file since 1952 was 1851, checked in by maronga, 8 years ago

last commit documented

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