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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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