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

Last change on this file since 1849 was 1849, checked in by hoffmann, 6 years ago

lpm_droplet_condensation improved, microphysics partially modularized

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