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

Last change on this file since 1691 was 1691, checked in by maronga, 6 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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