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

Last change on this file since 1734 was 1692, checked in by maronga, 9 years ago

last commit documented

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