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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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