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

Last change on this file since 2074 was 2032, checked in by knoop, 7 years ago

last commit documented

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