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

Last change on this file since 1671 was 1647, checked in by hoffmann, 9 years ago

last commit documented

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