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

Last change on this file since 2564 was 2522, checked in by schwenkel, 7 years ago

minor bugfix

  • Property svn:keywords set to Id
File size: 124.6 KB
Line 
1!> @file microphysics_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: microphysics_mod.f90 2522 2017-10-05 14:20:37Z Giersch $
27! Minor bugfix
28!
29! 2375 2017-08-29 14:10:28Z schwenkel
30! Improved aerosol initilization and some minor bugfixes
31! for droplet sedimenation
32!
33! 2318 2017-07-20 17:27:44Z suehring
34! Get topography top index via Function call
35!
36! 2317 2017-07-20 17:27:19Z suehring
37! s1 changed to log_sigma
38!
39! 2292 2017-06-20 09:51:42Z schwenkel
40! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
41! includes two more prognostic equations for cloud drop concentration (nc)
42! and cloud water content (qc).
43! - The process of activation is parameterized with a simple Twomey
44!   activion scheme or with considering solution and curvature
45!   effects (Khvorostyanov and Curry ,2006).
46! - The saturation adjustment scheme is replaced by the parameterization
47!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
48! - All other microphysical processes of Seifert and Beheng are used.
49!   Additionally, in those processes the reduction of cloud number concentration
50!   is considered.
51!
52! 2233 2017-05-30 18:08:54Z suehring
53!
54! 2232 2017-05-30 17:47:52Z suehring
55! Adjustments to new topography and surface concept
56!
57! 2155 2017-02-21 09:57:40Z hoffmann
58! Bugfix in the calculation of microphysical quantities on ghost points.
59!
60! 2031 2016-10-21 15:11:58Z knoop
61! renamed variable rho to rho_ocean
62!
63! 2000 2016-08-20 18:09:15Z knoop
64! Forced header and separation lines into 80 columns
65!
66! 1850 2016-04-08 13:29:27Z maronga
67! Module renamed
68! Adapted for modularization of microphysics.
69!
70! 1845 2016-04-08 08:29:13Z raasch
71! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
72! point (instead of one point above surface)
73!
74! 1831 2016-04-07 13:15:51Z hoffmann
75! turbulence renamed collision_turbulence,
76! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
77! avaialble for microphysics_kessler.
78!
79! 1822 2016-04-07 07:49:42Z hoffmann
80! Unused variables removed.
81! Kessler scheme integrated.
82!
83! 1691 2015-10-26 16:17:44Z maronga
84! Added new routine calc_precipitation_amount. The routine now allows to account
85! for precipitation due to sedimenation of cloud (fog) droplets
86!
87! 1682 2015-10-07 23:56:08Z knoop
88! Code annotations made doxygen readable
89!
90! 1646 2015-09-02 16:00:10Z hoffmann
91! Bugfix: Wrong computation of d_mean.
92!
93! 1361 2014-04-16 15:17:48Z hoffmann
94! Bugfix in sedimentation_rain: Index corrected.
95! Vectorized version of adjust_cloud added.
96! Little reformatting of the code.
97!
98! 1353 2014-04-08 15:21:23Z heinze
99! REAL constants provided with KIND-attribute
100!
101! 1346 2014-03-27 13:18:20Z heinze
102! Bugfix: REAL constants provided with KIND-attribute especially in call of
103! intrinsic function like MAX, MIN, SIGN
104!
105! 1334 2014-03-25 12:21:40Z heinze
106! Bugfix: REAL constants provided with KIND-attribute
107!
108! 1322 2014-03-20 16:38:49Z raasch
109! REAL constants defined as wp-kind
110!
111! 1320 2014-03-20 08:40:49Z raasch
112! ONLY-attribute added to USE-statements,
113! kind-parameters added to all INTEGER and REAL declaration statements,
114! kinds are defined in new module kinds,
115! comment fields (!:) to be used for variable explanations added to
116! all variable declaration statements
117!
118! 1241 2013-10-30 11:36:58Z heinze
119! hyp and rho_ocean have to be calculated at each time step if data from external
120! file LSF_DATA are used
121!
122! 1115 2013-03-26 18:16:16Z hoffmann
123! microphyical tendencies are calculated in microphysics_control in an optimized
124! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
125!
126! 1106 2013-03-04 05:31:38Z raasch
127! small changes in code formatting
128!
129! 1092 2013-02-02 11:24:22Z raasch
130! unused variables removed
131! file put under GPL
132!
133! 1065 2012-11-22 17:42:36Z hoffmann
134! Sedimentation process implemented according to Stevens and Seifert (2008).
135! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
136! and Stevens, 2010).
137!
138! 1053 2012-11-13 17:11:03Z hoffmann
139! initial revision
140!
141! Description:
142! ------------
143!> Calculate bilk cloud microphysics.
144!------------------------------------------------------------------------------!
145 MODULE microphysics_mod
146
147    USE  kinds
148
149    IMPLICIT NONE
150
151    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
152    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
153    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
154    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
155    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
156
157    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
158    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
159    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
160    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
161    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
162    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
163    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
164    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
165    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
166    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
167    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
168    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
169    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
170    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
171    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
172    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
173    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
174    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
175    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
176    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
177    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
178    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
179    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
180    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
181    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
182    REAL(wp) ::  kappa_rr = 60.7_wp        !< const. in collision kernel (kg-1/3)
183    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
184    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
185    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
186    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
187    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
188    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
189    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
190    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
191    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
192    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
193    REAL(wp) ::  xcmax = 2.6E-10_wp        !< maximum cloud drop size (kg) (~ 40µm)
194    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
195    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
196
197    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
198    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
199    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
200    REAL(wp) ::  dt_micro                        !< microphysics time step
201    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
202    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
203    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
204    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
205    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
206    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
207
208    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nc_1d  !<
209    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nr_1d  !<
210    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_1d  !<
211    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qc_1d  !<
212    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qr_1d  !<
213    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_1d   !<
214
215    SAVE
216
217    PRIVATE
218    PUBLIC microphysics_control, microphysics_init
219
220    PUBLIC cloud_water_sedimentation, collision_turbulence,                    &
221           curvature_solution_effects_bulk, c_sedimentation,                   &
222           dry_aerosol_radius, dt_precipitation,                               &
223           limiter_sedimentation, na_init, nc_const, sigma_gc, sigma_bulk,     &
224           ventilation_effect
225
226
227    INTERFACE microphysics_control
228       MODULE PROCEDURE microphysics_control
229       MODULE PROCEDURE microphysics_control_ij
230    END INTERFACE microphysics_control
231
232    INTERFACE adjust_cloud
233       MODULE PROCEDURE adjust_cloud
234       MODULE PROCEDURE adjust_cloud_ij
235    END INTERFACE adjust_cloud
236
237    INTERFACE activation
238       MODULE PROCEDURE activation
239       MODULE PROCEDURE activation_ij
240    END INTERFACE activation
241
242    INTERFACE condensation
243       MODULE PROCEDURE condensation
244       MODULE PROCEDURE condensation_ij
245    END INTERFACE condensation
246
247    INTERFACE autoconversion
248       MODULE PROCEDURE autoconversion
249       MODULE PROCEDURE autoconversion_ij
250    END INTERFACE autoconversion
251
252    INTERFACE autoconversion_kessler
253       MODULE PROCEDURE autoconversion_kessler
254       MODULE PROCEDURE autoconversion_kessler_ij
255    END INTERFACE autoconversion_kessler
256
257    INTERFACE accretion
258       MODULE PROCEDURE accretion
259       MODULE PROCEDURE accretion_ij
260    END INTERFACE accretion
261
262    INTERFACE selfcollection_breakup
263       MODULE PROCEDURE selfcollection_breakup
264       MODULE PROCEDURE selfcollection_breakup_ij
265    END INTERFACE selfcollection_breakup
266
267    INTERFACE evaporation_rain
268       MODULE PROCEDURE evaporation_rain
269       MODULE PROCEDURE evaporation_rain_ij
270    END INTERFACE evaporation_rain
271
272    INTERFACE sedimentation_cloud
273       MODULE PROCEDURE sedimentation_cloud
274       MODULE PROCEDURE sedimentation_cloud_ij
275    END INTERFACE sedimentation_cloud
276
277    INTERFACE sedimentation_rain
278       MODULE PROCEDURE sedimentation_rain
279       MODULE PROCEDURE sedimentation_rain_ij
280    END INTERFACE sedimentation_rain
281
282    INTERFACE calc_precipitation_amount
283       MODULE PROCEDURE calc_precipitation_amount
284       MODULE PROCEDURE calc_precipitation_amount_ij
285    END INTERFACE calc_precipitation_amount
286
287 CONTAINS
288!------------------------------------------------------------------------------!
289! Description:
290! ------------
291!> Initialization of bulk microphysics
292!------------------------------------------------------------------------------!
293    SUBROUTINE microphysics_init
294
295       USE arrays_3d,                                                          &
296           ONLY:  dzu
297
298       USE constants,                                                          &
299           ONLY:  pi
300
301       USE cloud_parameters,                                                   &
302           ONLY:  molecular_weight_of_solute, rho_l, rho_s, vanthoff
303
304       USE control_parameters,                                                 &
305           ONLY:  aerosol_nacl, aerosol_c3h4o4 , aerosol_nh4no3,               &
306                  microphysics_morrison, microphysics_seifert
307
308       USE indices,                                                            &
309           ONLY:  nzb, nzt
310
311       IMPLICIT NONE
312
313!
314!--    constant for the sedimentation of cloud water (2-moment cloud physics)
315       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
316                          )**( 2.0_wp / 3.0_wp ) *                             &
317                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
318
319!
320!--    Calculate timestep according to precipitation
321       IF ( microphysics_seifert )  THEN
322          dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
323                             w_precipitation
324       ENDIF
325
326!
327!--    Set constants for certain aerosol type
328       IF ( microphysics_morrison )  THEN
329          IF ( aerosol_nacl ) THEN
330             molecular_weight_of_solute = 0.05844_wp 
331             rho_s                      = 2165.0_wp
332             vanthoff                   = 2.0_wp
333          ELSEIF ( aerosol_c3h4o4 ) THEN
334             molecular_weight_of_solute = 0.10406_wp 
335             rho_s                      = 1600.0_wp
336             vanthoff                   = 1.37_wp
337          ELSEIF ( aerosol_nh4no3 ) THEN
338             molecular_weight_of_solute = 0.08004_wp 
339             rho_s                      = 1720.0_wp
340             vanthoff                   = 2.31_wp
341          ENDIF
342       ENDIF
343 
344!
345!--    Pre-calculate frequently calculated fractions of pi and rho_l
346       pirho_l  = pi * rho_l / 6.0_wp
347       dpirho_l = 1.0_wp / pirho_l
348
349!
350!--    Allocate 1D microphysics arrays
351       ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
352                  qc_1d(nzb:nzt+1) )
353
354       IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
355          ALLOCATE ( nc_1d(nzb:nzt+1) )
356       ENDIF
357
358       IF ( microphysics_seifert )  THEN
359          ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
360       ENDIF
361
362    END SUBROUTINE microphysics_init
363
364
365!------------------------------------------------------------------------------!
366! Description:
367! ------------
368!> Control of microphysics for all grid points
369!------------------------------------------------------------------------------!
370    SUBROUTINE microphysics_control
371
372       USE arrays_3d,                                                          &
373           ONLY:  hyp, pt_init, prr, zu
374
375       USE cloud_parameters,                                                   &
376           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
377
378       USE control_parameters,                                                 &
379           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
380                  intermediate_timestep_count, large_scale_forcing,            &
381                  lsf_surf, microphysics_kessler, microphysics_morrison,       &
382                  microphysics_seifert, pt_surface,                            &
383                  rho_surface,surface_pressure
384
385       USE indices,                                                            &
386           ONLY:  nzb, nzt
387
388       USE kinds
389
390       USE statistics,                                                         &
391           ONLY:  weight_pres
392
393       IMPLICIT NONE
394
395       INTEGER(iwp) ::  k                 !<
396
397       REAL(wp)     ::  t_surface         !<
398
399       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
400!
401!--       Calculate:
402!--       pt / t : ratio of potential and actual temperature (pt_d_t)
403!--       t / pt : ratio of actual and potential temperature (t_d_pt)
404!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
405          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
406          DO  k = nzb, nzt+1
407             hyp(k)    = surface_pressure * 100.0_wp * &
408                         ( ( t_surface - g / cp * zu(k) ) /                    &
409                         t_surface )**(1.0_wp / 0.286_wp)
410             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
411             t_d_pt(k) = 1.0_wp / pt_d_t(k)
412             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
413          ENDDO
414
415!
416!--       Compute reference density
417          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
418       ENDIF
419
420!
421!--    Compute length of time step
422       IF ( call_microphysics_at_all_substeps )  THEN
423          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
424       ELSE
425          dt_micro = dt_3d
426       ENDIF
427
428!
429!--    Reset precipitation rate
430       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
431
432!
433!--    Compute cloud physics
434       IF ( microphysics_kessler )  THEN
435
436          CALL autoconversion_kessler
437          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
438
439       ELSEIF ( microphysics_seifert )  THEN
440
441          CALL adjust_cloud
442          IF ( microphysics_morrison )  CALL activation
443          IF ( microphysics_morrison )  CALL condensation
444          CALL autoconversion
445          CALL accretion
446          CALL selfcollection_breakup
447          CALL evaporation_rain
448          CALL sedimentation_rain
449          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
450
451       ENDIF
452
453       CALL calc_precipitation_amount
454
455    END SUBROUTINE microphysics_control
456
457!------------------------------------------------------------------------------!
458! Description:
459! ------------
460!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
461!> evaporation of rain drops due to too small or too big weights
462!> of rain drops (Stevens and Seifert, 2008).
463!------------------------------------------------------------------------------!
464    SUBROUTINE adjust_cloud
465
466       USE arrays_3d,                                                          &
467           ONLY:  qc, nc, qr, nr
468
469       USE cloud_parameters,                                                   &
470           ONLY:  hyrho
471
472       USE control_parameters,                                                 &
473           ONLY:  microphysics_morrison
474
475       USE cpulog,                                                             &
476           ONLY:  cpu_log, log_point_s
477
478       USE indices,                                                            &
479           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
480
481       USE kinds
482
483       IMPLICIT NONE
484
485       INTEGER(iwp) ::  i                 !<
486       INTEGER(iwp) ::  j                 !<
487       INTEGER(iwp) ::  k                 !<
488
489       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
490
491       DO  i = nxlg, nxrg
492          DO  j = nysg, nyng
493             DO  k = nzb+1, nzt
494                IF ( qr(k,j,i) <= eps_sb )  THEN
495                   qr(k,j,i) = 0.0_wp
496                   nr(k,j,i) = 0.0_wp
497                ELSE
498                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
499                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
500                                       MERGE( 1.0_wp, 0.0_wp,                  &
501                                              BTEST( wall_flags_0(k,j,i), 0 ) )
502                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
503                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
504                                       MERGE( 1.0_wp, 0.0_wp,                  &
505                                              BTEST( wall_flags_0(k,j,i), 0 ) )
506                   ENDIF
507                ENDIF
508                IF ( microphysics_morrison ) THEN
509                   IF ( qc(k,j,i) <= eps_sb )  THEN
510                      qc(k,j,i) = 0.0_wp
511                      nc(k,j,i) = 0.0_wp
512                   ELSE
513                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
514                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
515                                          MERGE( 1.0_wp, 0.0_wp,               &
516                                              BTEST( wall_flags_0(k,j,i), 0 ) )
517                      ENDIF
518                   ENDIF
519                ENDIF
520             ENDDO
521          ENDDO
522       ENDDO
523
524       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
525
526    END SUBROUTINE adjust_cloud
527
528!------------------------------------------------------------------------------!
529! Description:
530! ------------
531!> Calculate number of activated condensation nucleii after simple activation
532!> scheme of Twomey, 1959.
533!------------------------------------------------------------------------------!
534    SUBROUTINE activation
535
536       USE arrays_3d,                                                          &
537           ONLY:  hyp, nc, nr, pt, q,  qc, qr
538
539       USE cloud_parameters,                                                   &
540           ONLY:  hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute,       &
541                  molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt,        &
542                  vanthoff
543
544       USE constants,                                                          &
545           ONLY:  pi
546
547       USE cpulog,                                                             &
548           ONLY:  cpu_log, log_point_s
549
550       USE indices,                                                            &
551           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
552
553       USE kinds
554
555       USE control_parameters,                                                 &
556          ONLY: simulated_time
557
558       IMPLICIT NONE
559
560       INTEGER(iwp) ::  i                 !<
561       INTEGER(iwp) ::  j                 !<
562       INTEGER(iwp) ::  k                 !<
563
564       REAL(wp)     ::  activ             !<
565       REAL(wp)     ::  afactor           !<
566       REAL(wp)     ::  alpha             !<
567       REAL(wp)     ::  beta_act          !<
568       REAL(wp)     ::  bfactor           !<
569       REAL(wp)     ::  e_s               !<
570       REAL(wp)     ::  k_act             !<
571       REAL(wp)     ::  n_act             !<
572       REAL(wp)     ::  n_ccn             !<
573       REAL(wp)     ::  q_s               !<
574       REAL(wp)     ::  s_0               !<
575       REAL(wp)     ::  sat               !<
576       REAL(wp)     ::  sat_max           !<
577       REAL(wp)     ::  sigma             !<
578       REAL(wp)     ::  sigma_act         !<
579       REAL(wp)     ::  t_int             !<
580       REAL(wp)     ::  t_l               !<
581
582       CALL cpu_log( log_point_s(65), 'activation', 'start' )
583
584       DO  i = nxlg, nxrg
585          DO  j = nysg, nyng
586             DO  k = nzb+1, nzt
587
588!
589!--             Actual liquid water temperature:
590                t_l = t_d_pt(k) * pt(k,j,i)
591
592!
593!--             Calculate actual temperature
594                t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
595!
596!--             Saturation vapor pressure at t_l:
597                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
598                                       ( t_l - 35.86_wp )                   &
599                                     )
600!
601!--             Computation of saturation humidity:
602                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
603                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
604                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
605                        ( 1.0_wp + alpha * q_s )
606
607!--             Supersaturation:
608                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
609
610!
611!--             Prescribe parameters for activation
612!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
613                k_act  = 0.7_wp
614                activ  = 0.0_wp
615
616
617                IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
618!
619!--                Compute the number of activated Aerosols
620!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
621                   n_act     = na_init * sat**k_act
622!
623!--                Compute the number of cloud droplets
624!--                (see: Morrison + Grabowski, 2007, JAS, 64)
625!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
626
627!
628!--                Compute activation rate after Khairoutdinov and Kogan
629!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
630                   sat_max = 1.0_wp / 100.0_wp
631                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
632                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
633                       dt_micro
634                ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
635!
636!--                Curvature effect (afactor) with surface tension
637!--                parameterization by Straka (2009)
638                   sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
639                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
640!
641!--                Solute effect (bfactor)
642                   bfactor = vanthoff * molecular_weight_of_water *            &
643                       rho_s / ( molecular_weight_of_solute * rho_l )
644
645!
646!--                Prescribe power index that describes the soluble fraction
647!--                of an aerosol particle (beta)
648!--                (see: Morrison + Grabowski, 2007, JAS, 64)
649                   beta_act  = 0.5_wp
650                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
651!
652!--                Calculate mean geometric supersaturation (s_0) with
653!--                parameterization by Khvorostyanov and Curry (2006)
654                   s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
655                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
656
657!
658!--                Calculate number of activated CCN as a function of
659!--                supersaturation and taking Koehler theory into account
660!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
661                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
662                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
663                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
664                ENDIF
665
666                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
667
668             ENDDO
669          ENDDO
670       ENDDO
671
672       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
673
674    END SUBROUTINE activation
675
676
677!------------------------------------------------------------------------------!
678! Description:
679! ------------
680!> Calculate condensation rate for cloud water content (after Khairoutdinov and
681!> Kogan, 2000).
682!------------------------------------------------------------------------------!
683    SUBROUTINE condensation
684
685       USE arrays_3d,                                                          &
686           ONLY:  hyp, nr, pt, q,  qc, qr, nc
687
688       USE cloud_parameters,                                                   &
689           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
690
691       USE constants,                                                          &
692           ONLY:  pi
693
694       USE cpulog,                                                             &
695           ONLY:  cpu_log, log_point_s
696
697       USE indices,                                                            &
698           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
699
700       USE kinds
701
702       USE control_parameters,                                                 &
703          ONLY: simulated_time
704
705       IMPLICIT NONE
706
707       INTEGER(iwp) ::  i                 !<
708       INTEGER(iwp) ::  j                 !<
709       INTEGER(iwp) ::  k                 !<
710
711       REAL(wp)     ::  alpha             !<
712       REAL(wp)     ::  cond              !<
713       REAL(wp)     ::  cond_max          !<
714       REAL(wp)     ::  dc                !<
715       REAL(wp)     ::  e_s               !<
716       REAL(wp)     ::  evap              !<
717       REAL(wp)     ::  evap_nc           !<
718       REAL(wp)     ::  g_fac             !<
719       REAL(wp)     ::  nc_0              !<
720       REAL(wp)     ::  q_s               !<
721       REAL(wp)     ::  sat               !<
722       REAL(wp)     ::  t_l               !<
723       REAL(wp)     ::  temp              !<
724       REAL(wp)     ::  xc                !<
725
726       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
727
728       DO  i = nxlg, nxrg
729          DO  j = nysg, nyng
730             DO  k = nzb+1, nzt
731!
732!--             Actual liquid water temperature:
733                t_l = t_d_pt(k) * pt(k,j,i)
734!
735!--             Saturation vapor pressure at t_l:
736                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /       &
737                                       ( t_l - 35.86_wp )                      &
738                                     )
739!
740!--             Computation of saturation humidity:
741                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
742                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
743                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /                  &
744                        ( 1.0_wp + alpha * q_s )
745
746!--             Supersaturation:
747                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
748
749!
750!--             Actual temperature:
751                temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
752
753                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
754                                    l_v / ( thermal_conductivity_l * temp )    &
755                                    + r_v * temp / ( diff_coeff_l * e_s )      &
756                                    )
757!
758!--             Mean weight of cloud drops
759                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
760                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
761!
762!--             Weight averaged diameter of cloud drops:
763                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
764!
765!--             Integral diameter of cloud drops
766                nc_0 = nc(k,j,i) * dc
767!
768!--             Condensation needs only to be calculated in supersaturated regions
769                IF ( sat > 0.0_wp )  THEN
770!
771!--                Condensation rate of cloud water content
772!--                after KK scheme.
773!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
774                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
775                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
776                   cond      = MIN( cond, cond_max / dt_micro )
777
778                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro
779                ELSEIF ( sat < 0.0_wp ) THEN
780                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
781                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
782
783                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro
784                ENDIF
785                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
786                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
787                ENDIF
788             ENDDO
789          ENDDO
790       ENDDO
791
792       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
793
794    END SUBROUTINE condensation
795
796
797!------------------------------------------------------------------------------!
798! Description:
799! ------------
800!> Autoconversion rate (Seifert and Beheng, 2006).
801!------------------------------------------------------------------------------!
802    SUBROUTINE autoconversion
803
804       USE arrays_3d,                                                          &
805           ONLY:  diss, dzu, nc, nr, qc, qr
806
807       USE cloud_parameters,                                                   &
808           ONLY:  hyrho
809
810       USE control_parameters,                                                 &
811           ONLY:  microphysics_morrison, rho_surface
812
813       USE cpulog,                                                             &
814           ONLY:  cpu_log, log_point_s
815
816       USE grid_variables,                                                     &
817           ONLY:  dx, dy
818
819       USE indices,                                                            &
820           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
821
822       USE kinds
823
824       IMPLICIT NONE
825
826       INTEGER(iwp) ::  i                 !<
827       INTEGER(iwp) ::  j                 !<
828       INTEGER(iwp) ::  k                 !<
829
830       REAL(wp)     ::  alpha_cc          !<
831       REAL(wp)     ::  autocon           !<
832       REAL(wp)     ::  dissipation       !<
833       REAL(wp)     ::  flag              !< flag to mask topography grid points
834       REAL(wp)     ::  k_au              !<
835       REAL(wp)     ::  l_mix             !<
836       REAL(wp)     ::  nc_auto           !<
837       REAL(wp)     ::  nu_c              !<
838       REAL(wp)     ::  phi_au            !<
839       REAL(wp)     ::  r_cc              !<
840       REAL(wp)     ::  rc                !<
841       REAL(wp)     ::  re_lambda         !<
842       REAL(wp)     ::  sigma_cc          !<
843       REAL(wp)     ::  tau_cloud         !<
844       REAL(wp)     ::  xc                !<
845
846       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
847
848       DO  i = nxlg, nxrg
849          DO  j = nysg, nyng
850             DO  k = nzb+1, nzt
851!
852!--             Predetermine flag to mask topography
853                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
854
855                nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
856
857                IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
858
859                   k_au = k_cc / ( 20.0_wp * x0 )
860!
861!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
862!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
863                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
864                                    qc(k,j,i) ), 0.0_wp )
865!
866!--                Universal function for autoconversion process
867!--                (Seifert and Beheng, 2006):
868                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
869                            ( 1.0_wp - tau_cloud**0.68_wp )**3
870!
871!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
872!--                (Use constant nu_c = 1.0_wp instead?)
873                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
874!
875!--                Mean weight of cloud droplets:
876                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
877!
878!--                Parameterized turbulence effects on autoconversion (Seifert,
879!--                Nuijens and Stevens, 2010)
880                   IF ( collision_turbulence )  THEN
881!
882!--                   Weight averaged radius of cloud droplets:
883                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
884
885                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
886                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
887                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
888!
889!--                   Mixing length (neglecting distance to ground and
890!--                   stratification)
891                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
892!
893!--                   Limit dissipation rate according to Seifert, Nuijens and
894!--                   Stevens (2010)
895                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
896!
897!--                   Compute Taylor-microscale Reynolds number:
898                      re_lambda = 6.0_wp / 11.0_wp *                           &
899                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
900                                  SQRT( 15.0_wp / kin_vis_air ) *              &
901                                  dissipation**( 1.0_wp / 6.0_wp )
902!
903!--                   The factor of 1.0E4 is needed to convert the dissipation
904!--                   rate from m2 s-3 to cm2 s-3.
905                      k_au = k_au * ( 1.0_wp +                                 &
906                                      dissipation * 1.0E4_wp *                 &
907                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
908                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
909                                                                      r_cc ) / &
910                                                        sigma_cc )**2          &
911                                                      ) + beta_cc              &
912                                      )                                        &
913                                    )
914                   ENDIF
915!
916!--                Autoconversion rate (Seifert and Beheng, 2006):
917                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
918                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
919                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
920                             rho_surface
921                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
922
923                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
924                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
925                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
926                                                              * flag
927                   IF ( microphysics_morrison ) THEN
928                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
929                                    autocon / x0 * hyrho(k) * dt_micro * flag )
930                   ENDIF
931
932                ENDIF
933
934             ENDDO
935          ENDDO
936       ENDDO
937
938       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
939
940    END SUBROUTINE autoconversion
941
942
943!------------------------------------------------------------------------------!
944! Description:
945! ------------
946!> Autoconversion process (Kessler, 1969).
947!------------------------------------------------------------------------------!
948    SUBROUTINE autoconversion_kessler
949
950       USE arrays_3d,                                                          &
951           ONLY:  dzw, pt, prr, q, qc
952
953       USE cloud_parameters,                                                   &
954           ONLY:  l_d_cp, pt_d_t
955
956       USE indices,                                                            &
957           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
958
959       USE kinds
960
961       USE surface_mod,                                                        &
962           ONLY:  get_topography_top_index
963
964
965       IMPLICIT NONE
966
967       INTEGER(iwp) ::  i      !<
968       INTEGER(iwp) ::  j      !<
969       INTEGER(iwp) ::  k      !<
970       INTEGER(iwp) ::  k_wall !< topgraphy top index
971
972       REAL(wp)    ::  dqdt_precip !<
973       REAL(wp)    ::  flag        !< flag to mask topography grid points
974
975       DO  i = nxlg, nxrg
976          DO  j = nysg, nyng
977!
978!--          Determine vertical index of topography top
979             k_wall = get_topography_top_index( j, i, 's' )
980             DO  k = nzb+1, nzt
981!
982!--             Predetermine flag to mask topography
983                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
984
985                IF ( qc(k,j,i) > ql_crit )  THEN
986                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
987                ELSE
988                   dqdt_precip = 0.0_wp
989                ENDIF
990
991                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
992                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
993                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
994                                        pt_d_t(k)              * flag
995
996!
997!--             Compute the rain rate (stored on surface grid point)
998                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
999
1000             ENDDO
1001          ENDDO
1002       ENDDO
1003
1004   END SUBROUTINE autoconversion_kessler
1005
1006
1007!------------------------------------------------------------------------------!
1008! Description:
1009! ------------
1010!> Accretion rate (Seifert and Beheng, 2006).
1011!------------------------------------------------------------------------------!
1012    SUBROUTINE accretion
1013
1014       USE arrays_3d,                                                          &
1015           ONLY:  diss, qc, qr, nc
1016
1017       USE cloud_parameters,                                                   &
1018           ONLY:  hyrho
1019
1020       USE control_parameters,                                                 &
1021           ONLY:  microphysics_morrison, rho_surface
1022
1023       USE cpulog,                                                             &
1024           ONLY:  cpu_log, log_point_s
1025
1026       USE indices,                                                            &
1027           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1028
1029       USE kinds
1030
1031       IMPLICIT NONE
1032
1033       INTEGER(iwp) ::  i                 !<
1034       INTEGER(iwp) ::  j                 !<
1035       INTEGER(iwp) ::  k                 !<
1036
1037       REAL(wp)     ::  accr              !<
1038       REAL(wp)     ::  flag              !< flag to mask topography grid points
1039       REAL(wp)     ::  k_cr              !<
1040       REAL(wp)     ::  nc_accr           !<
1041       REAL(wp)     ::  phi_ac            !<
1042       REAL(wp)     ::  tau_cloud         !<
1043       REAL(wp)     ::  xc                !<
1044
1045
1046       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
1047
1048       DO  i = nxlg, nxrg
1049          DO  j = nysg, nyng
1050             DO  k = nzb+1, nzt
1051!
1052!--             Predetermine flag to mask topography
1053                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1054
1055                nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
1056
1057                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
1058                                             .AND.  ( nc_accr > eps_mr ) ) THEN
1059!
1060!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
1061                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
1062!
1063!--                Universal function for accretion process (Seifert and
1064!--                Beheng, 2001):
1065                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1066
1067!
1068!--                Mean weight of cloud drops
1069                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
1070!
1071!--                Parameterized turbulence effects on autoconversion (Seifert,
1072!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1073!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1074                   IF ( collision_turbulence )  THEN
1075                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
1076                                       MIN( 600.0_wp,                          &
1077                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
1078                                     )
1079                   ELSE
1080                      k_cr = k_cr0
1081                   ENDIF
1082!
1083!--                Accretion rate (Seifert and Beheng, 2006):
1084                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
1085                          SQRT( rho_surface * hyrho(k) )
1086                   accr = MIN( accr, qc(k,j,i) / dt_micro )
1087
1088                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
1089                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
1090                   IF ( microphysics_morrison )  THEN
1091                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
1092                                         accr / xc * hyrho(k) * dt_micro * flag)
1093                   ENDIF
1094
1095                ENDIF
1096
1097             ENDDO
1098          ENDDO
1099       ENDDO
1100
1101       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
1102
1103    END SUBROUTINE accretion
1104
1105
1106!------------------------------------------------------------------------------!
1107! Description:
1108! ------------
1109!> Collisional breakup rate (Seifert, 2008).
1110!------------------------------------------------------------------------------!
1111    SUBROUTINE selfcollection_breakup
1112
1113       USE arrays_3d,                                                          &
1114           ONLY:  nr, qr
1115
1116       USE cloud_parameters,                                                   &
1117           ONLY:  hyrho
1118
1119       USE control_parameters,                                                 &
1120           ONLY:  rho_surface
1121
1122       USE cpulog,                                                             &
1123           ONLY:  cpu_log, log_point_s
1124
1125       USE indices,                                                            &
1126           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1127
1128       USE kinds
1129
1130       IMPLICIT NONE
1131
1132       INTEGER(iwp) ::  i                 !<
1133       INTEGER(iwp) ::  j                 !<
1134       INTEGER(iwp) ::  k                 !<
1135
1136       REAL(wp)     ::  breakup           !<
1137       REAL(wp)     ::  dr                !<
1138       REAL(wp)     ::  flag              !< flag to mask topography grid points
1139       REAL(wp)     ::  phi_br            !<
1140       REAL(wp)     ::  selfcoll          !<
1141
1142       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
1143
1144       DO  i = nxlg, nxrg
1145          DO  j = nysg, nyng
1146             DO  k = nzb+1, nzt
1147!
1148!--             Predetermine flag to mask topography
1149                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1150
1151                IF ( qr(k,j,i) > eps_sb )  THEN
1152!
1153!--                Selfcollection rate (Seifert and Beheng, 2001):
1154                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
1155                              SQRT( hyrho(k) * rho_surface )
1156!
1157!--                Weight averaged diameter of rain drops:
1158                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1159                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1160!
1161!--                Collisional breakup rate (Seifert, 2008):
1162                   IF ( dr >= 0.3E-3_wp )  THEN
1163                      phi_br  = k_br * ( dr - 1.1E-3_wp )
1164                      breakup = selfcoll * ( phi_br + 1.0_wp )
1165                   ELSE
1166                      breakup = 0.0_wp
1167                   ENDIF
1168
1169                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
1170                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
1171
1172                ENDIF
1173             ENDDO
1174          ENDDO
1175       ENDDO
1176
1177       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
1178
1179    END SUBROUTINE selfcollection_breakup
1180
1181
1182!------------------------------------------------------------------------------!
1183! Description:
1184! ------------
1185!> Evaporation of precipitable water. Condensation is neglected for
1186!> precipitable water.
1187!------------------------------------------------------------------------------!
1188    SUBROUTINE evaporation_rain
1189
1190       USE arrays_3d,                                                          &
1191           ONLY:  hyp, nr, pt, q,  qc, qr
1192
1193       USE cloud_parameters,                                                   &
1194           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
1195
1196       USE constants,                                                          &
1197           ONLY:  pi
1198
1199       USE cpulog,                                                             &
1200           ONLY:  cpu_log, log_point_s
1201
1202       USE indices,                                                            &
1203           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1204
1205       USE kinds
1206
1207       IMPLICIT NONE
1208
1209       INTEGER(iwp) ::  i                 !<
1210       INTEGER(iwp) ::  j                 !<
1211       INTEGER(iwp) ::  k                 !<
1212
1213       REAL(wp)     ::  alpha             !<
1214       REAL(wp)     ::  dr                !<
1215       REAL(wp)     ::  e_s               !<
1216       REAL(wp)     ::  evap              !<
1217       REAL(wp)     ::  evap_nr           !<
1218       REAL(wp)     ::  f_vent            !<
1219       REAL(wp)     ::  flag              !< flag to mask topography grid points
1220       REAL(wp)     ::  g_evap            !<
1221       REAL(wp)     ::  lambda_r          !<
1222       REAL(wp)     ::  mu_r              !<
1223       REAL(wp)     ::  mu_r_2            !<
1224       REAL(wp)     ::  mu_r_5d2          !<
1225       REAL(wp)     ::  nr_0              !<
1226       REAL(wp)     ::  q_s               !<
1227       REAL(wp)     ::  sat               !<
1228       REAL(wp)     ::  t_l               !<
1229       REAL(wp)     ::  temp              !<
1230       REAL(wp)     ::  xr                !<
1231
1232       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
1233
1234       DO  i = nxlg, nxrg
1235          DO  j = nysg, nyng
1236             DO  k = nzb+1, nzt
1237!
1238!--             Predetermine flag to mask topography
1239                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1240
1241                IF ( qr(k,j,i) > eps_sb )  THEN
1242!
1243!--                Actual liquid water temperature:
1244                   t_l = t_d_pt(k) * pt(k,j,i)
1245!
1246!--                Saturation vapor pressure at t_l:
1247                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
1248                                          ( t_l - 35.86_wp )                   &
1249                                        )
1250!
1251!--                Computation of saturation humidity:
1252                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1253                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1254                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
1255                           ( 1.0_wp + alpha * q_s )
1256!
1257!--                Supersaturation:
1258                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
1259!
1260!--                Evaporation needs only to be calculated in subsaturated regions
1261                   IF ( sat < 0.0_wp )  THEN
1262!
1263!--                   Actual temperature:
1264                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
1265
1266                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
1267                                          l_v / ( thermal_conductivity_l * temp ) &
1268                                          + r_v * temp / ( diff_coeff_l * e_s )   &
1269                                        )
1270!
1271!--                   Mean weight of rain drops
1272                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
1273!
1274!--                   Weight averaged diameter of rain drops:
1275                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1276!
1277!--                   Compute ventilation factor and intercept parameter
1278!--                   (Seifert and Beheng, 2006; Seifert, 2008):
1279                      IF ( ventilation_effect )  THEN
1280!
1281!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
1282!--                      2005; Stevens and Seifert, 2008):
1283                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
1284                                                          ( dr - 1.4E-3_wp ) ) )
1285!
1286!--                      Slope parameter of gamma distribution (Seifert, 2008):
1287                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
1288                                      ( mu_r + 1.0_wp )                        &
1289                                    )**( 1.0_wp / 3.0_wp ) / dr
1290
1291                         mu_r_2   = mu_r + 2.0_wp
1292                         mu_r_5d2 = mu_r + 2.5_wp
1293
1294                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
1295                                  lambda_r**( -mu_r_2 ) + b_vent *              &
1296                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
1297                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
1298                                  ( 1.0_wp -                                    &
1299                                    0.5_wp * ( b_term / a_term ) *              &
1300                                    ( lambda_r / ( c_term + lambda_r )          &
1301                                    )**mu_r_5d2 -                               &
1302                                    0.125_wp * ( b_term / a_term )**2 *         &
1303                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
1304                                    )**mu_r_5d2 -                               &
1305                                    0.0625_wp * ( b_term / a_term )**3 *        &
1306                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
1307                                    )**mu_r_5d2 -                               &
1308                                    0.0390625_wp * ( b_term / a_term )**4 *     &
1309                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
1310                                    )**mu_r_5d2                                 &
1311                                  )
1312
1313                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
1314                                  gamm( mu_r + 1.0_wp )
1315                      ELSE
1316                         f_vent = 1.0_wp
1317                         nr_0   = nr(k,j,i) * dr
1318                      ENDIF
1319   !
1320   !--                Evaporation rate of rain water content (Seifert and
1321   !--                Beheng, 2006):
1322                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
1323                                hyrho(k)
1324                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
1325                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
1326                                     -nr(k,j,i) / dt_micro )
1327
1328                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
1329                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
1330
1331                   ENDIF
1332                ENDIF
1333
1334             ENDDO
1335          ENDDO
1336       ENDDO
1337
1338       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
1339
1340    END SUBROUTINE evaporation_rain
1341
1342
1343!------------------------------------------------------------------------------!
1344! Description:
1345! ------------
1346!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1347!------------------------------------------------------------------------------!
1348    SUBROUTINE sedimentation_cloud
1349
1350       USE arrays_3d,                                                          &
1351           ONLY:  ddzu, dzu, nc, pt, prr, q, qc
1352
1353       USE cloud_parameters,                                                   &
1354           ONLY:  hyrho, l_d_cp, pt_d_t
1355
1356       USE control_parameters,                                                 &
1357           ONLY:  call_microphysics_at_all_substeps,                           &
1358                  intermediate_timestep_count, microphysics_morrison
1359
1360       USE cpulog,                                                             &
1361           ONLY:  cpu_log, log_point_s
1362
1363       USE indices,                                                            &
1364           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1365
1366       USE kinds
1367
1368       USE statistics,                                                         &
1369           ONLY:  weight_substep
1370
1371
1372       IMPLICIT NONE
1373
1374       INTEGER(iwp) ::  i             !<
1375       INTEGER(iwp) ::  j             !<
1376       INTEGER(iwp) ::  k             !<
1377
1378       REAL(wp) ::  flag              !< flag to mask topography grid points
1379       REAL(wp) ::  nc_sedi           !<
1380
1381       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
1382       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
1383
1384
1385       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
1386
1387       sed_qc(nzt+1) = 0.0_wp
1388       sed_nc(nzt+1) = 0.0_wp
1389
1390       DO  i = nxlg, nxrg
1391          DO  j = nysg, nyng
1392             DO  k = nzt, nzb+1, -1
1393!
1394!--             Predetermine flag to mask topography
1395                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1396                nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
1397
1398!
1399!--             Sedimentation fluxes for number concentration are only calculated
1400!--             for cloud_scheme = 'morrison'
1401                IF ( microphysics_morrison ) THEN
1402                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
1403                      sed_nc(k) = sed_qc_const *                               &
1404                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
1405                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
1406                   ELSE
1407                      sed_nc(k) = 0.0_wp
1408                   ENDIF
1409
1410                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
1411                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
1412                                  ) * flag
1413
1414                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
1415                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
1416                ENDIF
1417
1418                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
1419                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
1420                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
1421                                                                           flag
1422                ELSE
1423                   sed_qc(k) = 0.0_wp
1424                ENDIF
1425
1426                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
1427                                            dt_micro + sed_qc(k+1)             &
1428                               ) * flag
1429
1430                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
1431                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1432                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
1433                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1434                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
1435                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1436                                        pt_d_t(k) * dt_micro            * flag
1437
1438!
1439!--             Compute the precipitation rate due to cloud (fog) droplets
1440                IF ( call_microphysics_at_all_substeps )  THEN
1441                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
1442                                * weight_substep(intermediate_timestep_count)  &
1443                                * flag
1444                ELSE
1445                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
1446                ENDIF
1447
1448             ENDDO
1449          ENDDO
1450       ENDDO
1451
1452       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
1453
1454    END SUBROUTINE sedimentation_cloud
1455
1456
1457!------------------------------------------------------------------------------!
1458! Description:
1459! ------------
1460!> Computation of sedimentation flux. Implementation according to Stevens
1461!> and Seifert (2008). Code is based on UCLA-LES.
1462!------------------------------------------------------------------------------!
1463    SUBROUTINE sedimentation_rain
1464
1465       USE arrays_3d,                                                          &
1466           ONLY:  ddzu, dzu, nr, pt, prr, q, qr
1467
1468       USE cloud_parameters,                                                   &
1469           ONLY:  hyrho, l_d_cp, pt_d_t
1470
1471       USE control_parameters,                                                 &
1472           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
1473       USE cpulog,                                                             &
1474           ONLY:  cpu_log, log_point_s
1475
1476       USE indices,                                                            &
1477           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1478
1479       USE kinds
1480
1481       USE statistics,                                                         &
1482           ONLY:  weight_substep
1483
1484       USE surface_mod,                                                        &
1485           ONLY :  bc_h
1486
1487       IMPLICIT NONE
1488
1489       INTEGER(iwp) ::  i             !< running index x direction
1490       INTEGER(iwp) ::  j             !< running index y direction
1491       INTEGER(iwp) ::  k             !< running index z direction
1492       INTEGER(iwp) ::  k_run         !<
1493       INTEGER(iwp) ::  l             !< running index of surface type
1494       INTEGER(iwp) ::  m             !< running index surface elements
1495       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1496       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
1497
1498       REAL(wp)     ::  c_run                      !<
1499       REAL(wp)     ::  d_max                      !<
1500       REAL(wp)     ::  d_mean                     !<
1501       REAL(wp)     ::  d_min                      !<
1502       REAL(wp)     ::  dr                         !<
1503       REAL(wp)     ::  flux                       !<
1504       REAL(wp)     ::  flag                       !< flag to mask topography grid points
1505       REAL(wp)     ::  lambda_r                   !<
1506       REAL(wp)     ::  mu_r                       !<
1507       REAL(wp)     ::  z_run                      !<
1508
1509       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1510       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1511       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1512       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1513       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1514       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1515       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1516       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
1517
1518       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
1519
1520!
1521!--    Compute velocities
1522       DO  i = nxlg, nxrg
1523          DO  j = nysg, nyng
1524             DO  k = nzb+1, nzt
1525!
1526!--             Predetermine flag to mask topography
1527                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1528
1529                IF ( qr(k,j,i) > eps_sb )  THEN
1530!
1531!--                Weight averaged diameter of rain drops:
1532                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1533                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1534!
1535!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1536!--                Stevens and Seifert, 2008):
1537                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
1538                                                     ( dr - 1.4E-3_wp ) ) )
1539!
1540!--                Slope parameter of gamma distribution (Seifert, 2008):
1541                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1542                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1543
1544                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1545                                               a_term - b_term * ( 1.0_wp +    &
1546                                                  c_term /                     &
1547                                                  lambda_r )**( -1.0_wp *      &
1548                                                  ( mu_r + 1.0_wp ) )          &
1549                                              )                                &
1550                                ) * flag
1551
1552                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1553                                               a_term - b_term * ( 1.0_wp +    &
1554                                                  c_term /                     &
1555                                                  lambda_r )**( -1.0_wp *      &
1556                                                  ( mu_r + 4.0_wp ) )          &
1557                                             )                                 &
1558                                ) * flag
1559                ELSE
1560                   w_nr(k) = 0.0_wp
1561                   w_qr(k) = 0.0_wp
1562                ENDIF
1563             ENDDO
1564!
1565!--          Adjust boundary values using surface data type.
1566!--          Upward-facing
1567             surf_s = bc_h(0)%start_index(j,i)
1568             surf_e = bc_h(0)%end_index(j,i)
1569             DO  m = surf_s, surf_e
1570                k         = bc_h(0)%k(m)
1571                w_nr(k-1) = w_nr(k)
1572                w_qr(k-1) = w_qr(k)
1573             ENDDO
1574!
1575!--          Downward-facing
1576             surf_s = bc_h(1)%start_index(j,i)
1577             surf_e = bc_h(1)%end_index(j,i)
1578             DO  m = surf_s, surf_e
1579                k         = bc_h(1)%k(m)
1580                w_nr(k+1) = w_nr(k)
1581                w_qr(k+1) = w_qr(k)
1582             ENDDO
1583!
1584!--          Model top boundary value
1585             w_nr(nzt+1) = 0.0_wp
1586             w_qr(nzt+1) = 0.0_wp
1587!
1588!--          Compute Courant number
1589             DO  k = nzb+1, nzt
1590!
1591!--             Predetermine flag to mask topography
1592                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1593
1594                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
1595                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
1596                          dt_micro * ddzu(k) * flag
1597                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
1598                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
1599                          dt_micro * ddzu(k) * flag
1600             ENDDO
1601!
1602!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1603             IF ( limiter_sedimentation )  THEN
1604
1605                DO k = nzb+1, nzt
1606!
1607!--                Predetermine flag to mask topography
1608                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1609
1610                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
1611                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1612                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1613
1614                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1615                                                              2.0_wp * d_max,  &
1616                                                              ABS( d_mean ) )  &
1617                                                      * flag
1618
1619                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
1620                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1621                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1622
1623                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1624                                                              2.0_wp * d_max,  &
1625                                                              ABS( d_mean ) )
1626                ENDDO
1627
1628             ELSE
1629
1630                nr_slope = 0.0_wp
1631                qr_slope = 0.0_wp
1632
1633             ENDIF
1634
1635             sed_nr(nzt+1) = 0.0_wp
1636             sed_qr(nzt+1) = 0.0_wp
1637!
1638!--          Compute sedimentation flux
1639             DO  k = nzt, nzb+1, -1
1640!
1641!--             Predetermine flag to mask topography
1642                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1643!
1644!--             Sum up all rain drop number densities which contribute to the flux
1645!--             through k-1/2
1646                flux  = 0.0_wp
1647                z_run = 0.0_wp ! height above z(k)
1648                k_run = k
1649                c_run = MIN( 1.0_wp, c_nr(k) )
1650                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1651                   flux  = flux + hyrho(k_run) *                               &
1652                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
1653                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
1654                                              * flag
1655                   z_run = z_run + dzu(k_run) * flag
1656                   k_run = k_run + 1          * flag
1657                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
1658                                              * flag
1659                ENDDO
1660!
1661!--             It is not allowed to sediment more rain drop number density than
1662!--             available
1663                flux = MIN( flux,                                              &
1664                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
1665                            dt_micro                                           &
1666                          )
1667
1668                sed_nr(k) = flux / dt_micro * flag
1669                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
1670                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1671!
1672!--             Sum up all rain water content which contributes to the flux
1673!--             through k-1/2
1674                flux  = 0.0_wp
1675                z_run = 0.0_wp ! height above z(k)
1676                k_run = k
1677                c_run = MIN( 1.0_wp, c_qr(k) )
1678
1679                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1680
1681                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
1682                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
1683                                  0.5_wp ) * c_run * dzu(k_run) * flag
1684                   z_run = z_run + dzu(k_run)                   * flag
1685                   k_run = k_run + 1                            * flag
1686                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
1687                                                                * flag
1688
1689                ENDDO
1690!
1691!--             It is not allowed to sediment more rain water content than
1692!--             available
1693                flux = MIN( flux,                                              &
1694                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
1695                            dt_micro                                           &
1696                          )
1697
1698                sed_qr(k) = flux / dt_micro * flag
1699
1700                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
1701                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1702                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
1703                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1704                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
1705                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1706                                        pt_d_t(k) * dt_micro            * flag
1707!
1708!--             Compute the rain rate
1709                IF ( call_microphysics_at_all_substeps )  THEN
1710                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
1711                                * weight_substep(intermediate_timestep_count) &
1712                                * flag
1713                ELSE
1714                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
1715                ENDIF
1716
1717             ENDDO
1718          ENDDO
1719       ENDDO
1720
1721       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1722
1723    END SUBROUTINE sedimentation_rain
1724
1725
1726!------------------------------------------------------------------------------!
1727! Description:
1728! ------------
1729!> Computation of the precipitation amount due to gravitational settling of
1730!> rain and cloud (fog) droplets
1731!------------------------------------------------------------------------------!
1732    SUBROUTINE calc_precipitation_amount
1733
1734       USE arrays_3d,                                                          &
1735           ONLY:  precipitation_amount, prr
1736
1737       USE cloud_parameters,                                                   &
1738           ONLY:  hyrho
1739
1740       USE control_parameters,                                                 &
1741           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
1742                  intermediate_timestep_count, intermediate_timestep_count_max,&
1743                  precipitation_amount_interval, time_do2d_xy
1744
1745       USE indices,                                                            &
1746           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
1747
1748       USE kinds
1749
1750       USE surface_mod,                                                        &
1751           ONLY :  bc_h
1752
1753       IMPLICIT NONE
1754
1755       INTEGER(iwp) ::  i             !< running index x direction
1756       INTEGER(iwp) ::  j             !< running index y direction
1757       INTEGER(iwp) ::  k             !< running index y direction
1758       INTEGER(iwp) ::  m             !< running index surface elements
1759       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1760       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
1761
1762       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1763            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
1764            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1765       THEN
1766!
1767!--       Run over all upward-facing surface elements, i.e. non-natural,
1768!--       natural and urban
1769          DO  m = 1, bc_h(0)%ns
1770             i = bc_h(0)%i(m)
1771             j = bc_h(0)%j(m)
1772             k = bc_h(0)%k(m)
1773             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
1774                                               prr(k,j,i) * hyrho(k) * dt_3d
1775          ENDDO
1776
1777       ENDIF
1778
1779    END SUBROUTINE calc_precipitation_amount
1780
1781
1782!------------------------------------------------------------------------------!
1783! Description:
1784! ------------
1785!> Control of microphysics for grid points i,j
1786!------------------------------------------------------------------------------!
1787
1788    SUBROUTINE microphysics_control_ij( i, j )
1789
1790       USE arrays_3d,                                                          &
1791           ONLY:  hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
1792
1793       USE cloud_parameters,                                                   &
1794           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
1795
1796       USE control_parameters,                                                 &
1797           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
1798                  intermediate_timestep_count, large_scale_forcing,            &
1799                  lsf_surf, microphysics_morrison, microphysics_seifert,       &
1800                  microphysics_kessler, pt_surface, rho_surface,               &
1801                  surface_pressure
1802
1803       USE indices,                                                            &
1804           ONLY:  nzb, nzt
1805
1806       USE kinds
1807
1808       USE statistics,                                                         &
1809           ONLY:  weight_pres
1810
1811       IMPLICIT NONE
1812
1813       INTEGER(iwp) ::  i                 !<
1814       INTEGER(iwp) ::  j                 !<
1815       INTEGER(iwp) ::  k                 !<
1816
1817       REAL(wp)     ::  t_surface         !<
1818
1819       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
1820!
1821!--       Calculate:
1822!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1823!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1824!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1825          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1826          DO  k = nzb, nzt+1
1827             hyp(k)    = surface_pressure * 100.0_wp * &
1828                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1829             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1830             t_d_pt(k) = 1.0_wp / pt_d_t(k)
1831             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
1832          ENDDO
1833!
1834!--       Compute reference density
1835          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1836       ENDIF
1837
1838!
1839!--    Compute length of time step
1840       IF ( call_microphysics_at_all_substeps )  THEN
1841          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1842       ELSE
1843          dt_micro = dt_3d
1844       ENDIF
1845
1846!
1847!--    Use 1d arrays
1848       q_1d(:)  = q(:,j,i)
1849       pt_1d(:) = pt(:,j,i)
1850       qc_1d(:) = qc(:,j,i)
1851       IF ( microphysics_seifert )  THEN
1852          qr_1d(:) = qr(:,j,i)
1853          nr_1d(:) = nr(:,j,i)
1854       ENDIF
1855       IF ( microphysics_morrison )  THEN
1856          nc_1d(:) = nc(:,j,i)
1857       ENDIF
1858
1859
1860!
1861!--    Reset precipitation rate
1862       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1863
1864!
1865!--    Compute cloud physics
1866       IF( microphysics_kessler )  THEN
1867
1868          CALL autoconversion_kessler( i,j )
1869          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1870
1871       ELSEIF ( microphysics_seifert )  THEN
1872
1873          CALL adjust_cloud( i,j )
1874          IF ( microphysics_morrison ) CALL activation( i,j )
1875          IF ( microphysics_morrison ) CALL condensation( i,j )
1876          CALL autoconversion( i,j )
1877          CALL accretion( i,j )
1878          CALL selfcollection_breakup( i,j )
1879          CALL evaporation_rain( i,j )
1880          CALL sedimentation_rain( i,j )
1881          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1882
1883       ENDIF
1884
1885       CALL calc_precipitation_amount( i,j )
1886
1887!
1888!--    Store results on the 3d arrays
1889       q(:,j,i)  = q_1d(:)
1890       pt(:,j,i) = pt_1d(:)
1891       IF ( microphysics_morrison ) THEN
1892          qc(:,j,i) = qc_1d(:)
1893          nc(:,j,i) = nc_1d(:)
1894       ENDIF
1895       IF ( microphysics_seifert )  THEN
1896          qr(:,j,i) = qr_1d(:)
1897          nr(:,j,i) = nr_1d(:)
1898       ENDIF
1899
1900    END SUBROUTINE microphysics_control_ij
1901
1902!------------------------------------------------------------------------------!
1903! Description:
1904! ------------
1905!> Adjust number of raindrops to avoid nonlinear effects in
1906!> sedimentation and evaporation of rain drops due to too small or
1907!> too big weights of rain drops (Stevens and Seifert, 2008).
1908!> The same procedure is applied to cloud droplets if they are determined
1909!> prognostically. Call for grid point i,j
1910!------------------------------------------------------------------------------!
1911    SUBROUTINE adjust_cloud_ij( i, j )
1912
1913       USE cloud_parameters,                                                   &
1914           ONLY:  hyrho
1915
1916       USE control_parameters,                                                 &
1917           ONLY: microphysics_morrison
1918
1919       USE indices,                                                            &
1920           ONLY:  nzb, nzt, wall_flags_0
1921
1922       USE kinds
1923
1924       IMPLICIT NONE
1925
1926       INTEGER(iwp) ::  i                 !<
1927       INTEGER(iwp) ::  j                 !<
1928       INTEGER(iwp) ::  k                 !<
1929
1930       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
1931
1932       DO  k = nzb+1, nzt
1933!
1934!--       Predetermine flag to mask topography
1935          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1936
1937          IF ( qr_1d(k) <= eps_sb )  THEN
1938             qr_1d(k) = 0.0_wp
1939             nr_1d(k) = 0.0_wp
1940          ELSE
1941!
1942!--          Adjust number of raindrops to avoid nonlinear effects in
1943!--          sedimentation and evaporation of rain drops due to too small or
1944!--          too big weights of rain drops (Stevens and Seifert, 2008).
1945             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
1946                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
1947             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
1948                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
1949             ENDIF
1950
1951          ENDIF
1952
1953          IF ( microphysics_morrison ) THEN
1954             IF ( qc_1d(k) <= eps_sb )  THEN
1955                qc_1d(k) = 0.0_wp
1956                nc_1d(k) = 0.0_wp
1957             ELSE
1958                IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) )  THEN
1959                   nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
1960                ENDIF
1961             ENDIF
1962          ENDIF
1963
1964       ENDDO
1965
1966    END SUBROUTINE adjust_cloud_ij
1967
1968!------------------------------------------------------------------------------!
1969! Description:
1970! ------------
1971!> Calculate number of activated condensation nucleii after simple activation
1972!> scheme of Twomey, 1959.
1973!------------------------------------------------------------------------------!
1974    SUBROUTINE activation_ij( i, j )
1975
1976       USE arrays_3d,                                                          &
1977           ONLY:  hyp, nr, pt, q,  qc, qr, nc
1978
1979       USE cloud_parameters,                                                   &
1980           ONLY:  hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute,       &
1981                  molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt,        &
1982                  vanthoff
1983
1984       USE constants,                                                          &
1985           ONLY:  pi
1986
1987       USE cpulog,                                                             &
1988           ONLY:  cpu_log, log_point_s
1989
1990       USE indices,                                                            &
1991           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
1992
1993       USE kinds
1994
1995       USE control_parameters,                                                 &
1996          ONLY: simulated_time
1997
1998       IMPLICIT NONE
1999
2000       INTEGER(iwp) ::  i                 !<
2001       INTEGER(iwp) ::  j                 !<
2002       INTEGER(iwp) ::  k                 !<
2003
2004       REAL(wp)     ::  activ             !<
2005       REAL(wp)     ::  afactor           !<
2006       REAL(wp)     ::  alpha             !<
2007       REAL(wp)     ::  beta_act          !<
2008       REAL(wp)     ::  bfactor           !<
2009       REAL(wp)     ::  e_s               !<
2010       REAL(wp)     ::  k_act             !<
2011       REAL(wp)     ::  n_act             !<
2012       REAL(wp)     ::  n_ccn             !<
2013       REAL(wp)     ::  q_s               !<
2014       REAL(wp)     ::  s_0               !<
2015       REAL(wp)     ::  sat               !<
2016       REAL(wp)     ::  sat_max           !<
2017       REAL(wp)     ::  sigma             !<
2018       REAL(wp)     ::  sigma_act         !<
2019       REAL(wp)     ::  t_int             !<
2020       REAL(wp)     ::  t_l               !<
2021
2022       DO  k = nzb+1, nzt
2023!
2024!--       Actual liquid water temperature:
2025          t_l = t_d_pt(k) * pt_1d(k)
2026
2027!
2028!--       Calculate actual temperature
2029          t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
2030!
2031!--             Saturation vapor pressure at t_l:
2032          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
2033                                 ( t_l - 35.86_wp )                            &
2034                                 )
2035!
2036!--       Computation of saturation humidity:
2037          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2038          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2039          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /                         &
2040                        ( 1.0_wp + alpha * q_s )
2041
2042!--       Supersaturation:
2043          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2044
2045!
2046!--       Prescribe parameters for activation
2047!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
2048          k_act  = 0.7_wp
2049          activ  = 0.0_wp
2050
2051          IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
2052!
2053!--          Compute the number of activated Aerosols
2054!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
2055             n_act     = na_init * sat**k_act
2056!
2057!--          Compute the number of cloud droplets
2058!--          (see: Morrison + Grabowski, 2007, JAS, 64)
2059!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
2060
2061!
2062!--          Compute activation rate after Khairoutdinov and Kogan
2063!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2064             sat_max = 0.8_wp / 100.0_wp
2065             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              &
2066                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          &
2067                  dt_micro
2068
2069             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
2070          ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
2071!
2072!--          Curvature effect (afactor) with surface tension
2073!--          parameterization by Straka (2009)
2074             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
2075             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
2076!
2077!--          Solute effect (bfactor)
2078             bfactor = vanthoff * molecular_weight_of_water *                  &
2079                  rho_s / ( molecular_weight_of_solute * rho_l )
2080
2081!
2082!--          Prescribe power index that describes the soluble fraction
2083!--          of an aerosol particle (beta).
2084!--          (see: Morrison + Grabowski, 2007, JAS, 64)
2085             beta_act  = 0.5_wp
2086             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2087!
2088!--          Calculate mean geometric supersaturation (s_0) with
2089!--          parameterization by Khvorostyanov and Curry (2006)
2090             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
2091               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
2092
2093!
2094!--          Calculate number of activated CCN as a function of
2095!--          supersaturation and taking Koehler theory into account
2096!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2097             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
2098                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2099             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
2100
2101             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
2102          ENDIF
2103
2104       ENDDO
2105
2106    END SUBROUTINE activation_ij
2107
2108!------------------------------------------------------------------------------!
2109! Description:
2110! ------------
2111!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2112!> Kogan, 2000).
2113!------------------------------------------------------------------------------!
2114    SUBROUTINE condensation_ij( i, j )
2115
2116       USE arrays_3d,                                                          &
2117           ONLY:  hyp, nr, pt, q,  qc, qr, nc
2118
2119       USE cloud_parameters,                                                   &
2120           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
2121
2122       USE constants,                                                          &
2123           ONLY:  pi
2124
2125       USE cpulog,                                                             &
2126           ONLY:  cpu_log, log_point_s
2127
2128       USE indices,                                                            &
2129           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
2130
2131       USE kinds
2132
2133       USE control_parameters,                                                 &
2134          ONLY: simulated_time
2135
2136       IMPLICIT NONE
2137
2138       INTEGER(iwp) ::  i                 !<
2139       INTEGER(iwp) ::  j                 !<
2140       INTEGER(iwp) ::  k                 !<
2141
2142       REAL(wp)     ::  alpha             !<
2143       REAL(wp)     ::  cond              !<
2144       REAL(wp)     ::  cond_max          !<
2145       REAL(wp)     ::  dc                !<
2146       REAL(wp)     ::  e_s               !<
2147       REAL(wp)     ::  evap              !<
2148       REAL(wp)     ::  evap_nc           !<
2149       REAL(wp)     ::  g_fac             !<
2150       REAL(wp)     ::  nc_0              !<
2151       REAL(wp)     ::  q_s               !<
2152       REAL(wp)     ::  sat               !<
2153       REAL(wp)     ::  t_l               !<
2154       REAL(wp)     ::  temp              !<
2155       REAL(wp)     ::  xc                !<
2156
2157
2158       DO  k = nzb+1, nzt
2159!
2160!--       Actual liquid water temperature:
2161          t_l = t_d_pt(k) * pt_1d(k)
2162!
2163!--             Saturation vapor pressure at t_l:
2164          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
2165                                 ( t_l - 35.86_wp )                   &
2166                                 )
2167!
2168!--       Computation of saturation humidity:
2169          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2170          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2171          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /               &
2172                        ( 1.0_wp + alpha * q_s )
2173
2174!--       Supersaturation:
2175          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2176
2177
2178!
2179!--       Actual temperature:
2180          temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
2181
2182          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
2183                              l_v / ( thermal_conductivity_l * temp )    &
2184                              + r_v * temp / ( diff_coeff_l * e_s )      &
2185                            )
2186!
2187!--       Mean weight of cloud drops
2188          IF ( nc_1d(k) <= 0.0_wp) CYCLE
2189          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
2190!
2191!--       Weight averaged diameter of cloud drops:
2192          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2193!
2194!--       Integral diameter of cloud drops
2195          nc_0 = nc_1d(k) * dc
2196!
2197!--       Condensation needs only to be calculated in supersaturated regions
2198          IF ( sat > 0.0_wp )  THEN
2199!
2200!--          Condensation rate of cloud water content
2201!--          after KK scheme.
2202!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2203             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2204             cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
2205             cond      = MIN( cond, cond_max / dt_micro )
2206
2207             qc_1d(k) = qc_1d(k) + cond * dt_micro
2208          ELSEIF ( sat < 0.0_wp ) THEN
2209             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2210             evap      = MAX( evap, -qc_1d(k) / dt_micro )
2211
2212             qc_1d(k) = qc_1d(k) + evap * dt_micro
2213          ENDIF
2214       ENDDO
2215
2216    END SUBROUTINE condensation_ij
2217
2218
2219!------------------------------------------------------------------------------!
2220! Description:
2221! ------------
2222!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
2223!------------------------------------------------------------------------------!
2224    SUBROUTINE autoconversion_ij( i, j )
2225
2226       USE arrays_3d,                                                          &
2227           ONLY:  diss, dzu
2228
2229       USE cloud_parameters,                                                   &
2230           ONLY:  hyrho
2231
2232       USE control_parameters,                                                 &
2233           ONLY:  microphysics_morrison, rho_surface
2234
2235       USE grid_variables,                                                     &
2236           ONLY:  dx, dy
2237
2238       USE indices,                                                            &
2239           ONLY:  nzb, nzt, wall_flags_0
2240
2241       USE kinds
2242
2243       IMPLICIT NONE
2244
2245       INTEGER(iwp) ::  i                 !<
2246       INTEGER(iwp) ::  j                 !<
2247       INTEGER(iwp) ::  k                 !<
2248
2249       REAL(wp)     ::  alpha_cc          !<
2250       REAL(wp)     ::  autocon           !<
2251       REAL(wp)     ::  dissipation       !<
2252       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
2253       REAL(wp)     ::  k_au              !<
2254       REAL(wp)     ::  l_mix             !<
2255       REAL(wp)     ::  nc_auto           !<
2256       REAL(wp)     ::  nu_c              !<
2257       REAL(wp)     ::  phi_au            !<
2258       REAL(wp)     ::  r_cc              !<
2259       REAL(wp)     ::  rc                !<
2260       REAL(wp)     ::  re_lambda         !<
2261       REAL(wp)     ::  sigma_cc          !<
2262       REAL(wp)     ::  tau_cloud         !<
2263       REAL(wp)     ::  xc                !<
2264
2265       DO  k = nzb+1, nzt
2266!
2267!--       Predetermine flag to mask topography
2268          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2269          nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
2270
2271          IF ( qc_1d(k) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
2272
2273             k_au = k_cc / ( 20.0_wp * x0 )
2274!
2275!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
2276!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
2277             tau_cloud = MAX( 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ),     & 
2278                              0.0_wp )
2279!
2280!--          Universal function for autoconversion process
2281!--          (Seifert and Beheng, 2006):
2282             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
2283!
2284!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
2285!--          (Use constant nu_c = 1.0_wp instead?)
2286             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
2287!
2288!--          Mean weight of cloud droplets:
2289             xc = hyrho(k) * qc_1d(k) / nc_auto
2290!
2291!--          Parameterized turbulence effects on autoconversion (Seifert,
2292!--          Nuijens and Stevens, 2010)
2293             IF ( collision_turbulence )  THEN
2294!
2295!--             Weight averaged radius of cloud droplets:
2296                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2297
2298                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2299                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2300                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2301!
2302!--             Mixing length (neglecting distance to ground and stratification)
2303                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2304!
2305!--             Limit dissipation rate according to Seifert, Nuijens and
2306!--             Stevens (2010)
2307                dissipation = MIN( 0.06_wp, diss(k,j,i) )
2308!
2309!--             Compute Taylor-microscale Reynolds number:
2310                re_lambda = 6.0_wp / 11.0_wp *                                 &
2311                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
2312                            SQRT( 15.0_wp / kin_vis_air ) *                    &
2313                            dissipation**( 1.0_wp / 6.0_wp )
2314!
2315!--             The factor of 1.0E4 is needed to convert the dissipation rate
2316!--             from m2 s-3 to cm2 s-3.
2317                k_au = k_au * ( 1.0_wp +                                       &
2318                                dissipation * 1.0E4_wp *                       &
2319                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
2320                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
2321                                                  sigma_cc )**2                &
2322                                                ) + beta_cc                    &
2323                                )                                              &
2324                              )
2325             ENDIF
2326!
2327!--          Autoconversion rate (Seifert and Beheng, 2006):
2328             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
2329                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
2330                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
2331                       rho_surface
2332             autocon = MIN( autocon, qc_1d(k) / dt_micro )
2333
2334             qr_1d(k) = qr_1d(k) + autocon * dt_micro                 * flag
2335             qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
2336             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
2337             IF ( microphysics_morrison ) THEN
2338                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                &
2339                                    autocon / x0 * hyrho(k) * dt_micro * flag )
2340             ENDIF
2341
2342          ENDIF
2343
2344       ENDDO
2345
2346    END SUBROUTINE autoconversion_ij
2347
2348!------------------------------------------------------------------------------!
2349! Description:
2350! ------------
2351!> Autoconversion process (Kessler, 1969).
2352!------------------------------------------------------------------------------!
2353    SUBROUTINE autoconversion_kessler_ij( i, j )
2354
2355       USE arrays_3d,                                                          &
2356           ONLY:  dzw, prr
2357
2358       USE cloud_parameters,                                                   &
2359           ONLY:  l_d_cp, pt_d_t
2360
2361       USE indices,                                                            &
2362           ONLY:  nzb, nzt, wall_flags_0
2363
2364       USE kinds
2365
2366       USE surface_mod,                                                        &
2367           ONLY:  get_topography_top_index
2368
2369
2370       IMPLICIT NONE
2371
2372       INTEGER(iwp) ::  i      !<
2373       INTEGER(iwp) ::  j      !<
2374       INTEGER(iwp) ::  k      !<
2375       INTEGER(iwp) ::  k_wall !< topography top index
2376
2377       REAL(wp)    ::  dqdt_precip !<
2378       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
2379
2380!
2381!--    Determine vertical index of topography top
2382       k_wall = get_topography_top_index( j, i, 's' )
2383       DO  k = nzb+1, nzt
2384!
2385!--       Predetermine flag to mask topography
2386          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2387
2388          IF ( qc_1d(k) > ql_crit )  THEN
2389             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
2390          ELSE
2391             dqdt_precip = 0.0_wp
2392          ENDIF
2393
2394          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
2395          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro * flag
2396          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
2397
2398!
2399!--       Compute the rain rate (stored on surface grid point)
2400          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2401
2402       ENDDO
2403
2404    END SUBROUTINE autoconversion_kessler_ij
2405
2406!------------------------------------------------------------------------------!
2407! Description:
2408! ------------
2409!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
2410!------------------------------------------------------------------------------!
2411    SUBROUTINE accretion_ij( i, j )
2412
2413       USE arrays_3d,                                                          &
2414           ONLY:  diss
2415
2416       USE cloud_parameters,                                                   &
2417           ONLY:  hyrho
2418
2419       USE control_parameters,                                                 &
2420           ONLY:  microphysics_morrison, rho_surface
2421
2422       USE indices,                                                            &
2423           ONLY:  nzb, nzt, wall_flags_0
2424
2425       USE kinds
2426
2427       IMPLICIT NONE
2428
2429       INTEGER(iwp) ::  i                 !<
2430       INTEGER(iwp) ::  j                 !<
2431       INTEGER(iwp) ::  k                 !<
2432
2433       REAL(wp)     ::  accr              !<
2434       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
2435       REAL(wp)     ::  k_cr              !<
2436       REAL(wp)     ::  nc_accr           !<
2437       REAL(wp)     ::  phi_ac            !<
2438       REAL(wp)     ::  tau_cloud         !<
2439       REAL(wp)     ::  xc                !<
2440
2441
2442       DO  k = nzb+1, nzt
2443!
2444!--       Predetermine flag to mask topography
2445          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2446          nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
2447
2448          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb )  .AND.  &
2449               ( nc_accr > eps_mr ) )  THEN
2450!
2451!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
2452             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
2453!
2454!--          Universal function for accretion process
2455!--          (Seifert and Beheng, 2001):
2456             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2457
2458!
2459!--          Mean weight of cloud drops
2460             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
2461!
2462!--          Parameterized turbulence effects on autoconversion (Seifert,
2463!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2464!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2465             IF ( collision_turbulence )  THEN
2466                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
2467                                 MIN( 600.0_wp,                                &
2468                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
2469                               )
2470             ELSE
2471                k_cr = k_cr0
2472             ENDIF
2473!
2474!--          Accretion rate (Seifert and Beheng, 2006):
2475             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                      &
2476                    SQRT( rho_surface * hyrho(k) )
2477             accr = MIN( accr, qc_1d(k) / dt_micro )
2478
2479             qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
2480             qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
2481             IF ( microphysics_morrison )  THEN
2482                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc *               &
2483                                             hyrho(k) * dt_micro * flag        &
2484                                         )
2485             ENDIF
2486
2487
2488          ENDIF
2489
2490       ENDDO
2491
2492    END SUBROUTINE accretion_ij
2493
2494
2495!------------------------------------------------------------------------------!
2496! Description:
2497! ------------
2498!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
2499!------------------------------------------------------------------------------!
2500    SUBROUTINE selfcollection_breakup_ij( i, j )
2501
2502       USE cloud_parameters,                                                   &
2503           ONLY:  hyrho
2504
2505       USE control_parameters,                                                 &
2506           ONLY:  rho_surface
2507
2508       USE indices,                                                            &
2509           ONLY:  nzb, nzt, wall_flags_0
2510
2511       USE kinds
2512
2513       IMPLICIT NONE
2514
2515       INTEGER(iwp) ::  i                 !<
2516       INTEGER(iwp) ::  j                 !<
2517       INTEGER(iwp) ::  k                 !<
2518
2519       REAL(wp)     ::  breakup           !<
2520       REAL(wp)     ::  dr                !<
2521       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
2522       REAL(wp)     ::  phi_br            !<
2523       REAL(wp)     ::  selfcoll          !<
2524
2525       DO  k = nzb+1, nzt
2526!
2527!--       Predetermine flag to mask topography
2528          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2529
2530          IF ( qr_1d(k) > eps_sb )  THEN
2531!
2532!--          Selfcollection rate (Seifert and Beheng, 2001):
2533             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
2534!
2535!--          Weight averaged diameter of rain drops:
2536             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
2537!
2538!--          Collisional breakup rate (Seifert, 2008):
2539             IF ( dr >= 0.3E-3_wp )  THEN
2540                phi_br  = k_br * ( dr - 1.1E-3_wp )
2541                breakup = selfcoll * ( phi_br + 1.0_wp )
2542             ELSE
2543                breakup = 0.0_wp
2544             ENDIF
2545
2546             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
2547             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
2548
2549          ENDIF
2550       ENDDO
2551
2552    END SUBROUTINE selfcollection_breakup_ij
2553
2554
2555!------------------------------------------------------------------------------!
2556! Description:
2557! ------------
2558!> Evaporation of precipitable water. Condensation is neglected for
2559!> precipitable water. Call for grid point i,j
2560!------------------------------------------------------------------------------!
2561    SUBROUTINE evaporation_rain_ij( i, j )
2562
2563       USE arrays_3d,                                                          &
2564           ONLY:  hyp
2565
2566       USE cloud_parameters,                                                   &
2567           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
2568
2569       USE constants,                                                          &
2570           ONLY:  pi
2571
2572       USE indices,                                                            &
2573           ONLY:  nzb, nzt, wall_flags_0
2574
2575       USE kinds
2576
2577       IMPLICIT NONE
2578
2579       INTEGER(iwp) ::  i                 !<
2580       INTEGER(iwp) ::  j                 !<
2581       INTEGER(iwp) ::  k                 !<
2582
2583       REAL(wp)     ::  alpha             !<
2584       REAL(wp)     ::  dr                !<
2585       REAL(wp)     ::  e_s               !<
2586       REAL(wp)     ::  evap              !<
2587       REAL(wp)     ::  evap_nr           !<
2588       REAL(wp)     ::  f_vent            !<
2589       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
2590       REAL(wp)     ::  g_evap            !<
2591       REAL(wp)     ::  lambda_r          !<
2592       REAL(wp)     ::  mu_r              !<
2593       REAL(wp)     ::  mu_r_2            !<
2594       REAL(wp)     ::  mu_r_5d2          !<
2595       REAL(wp)     ::  nr_0              !<
2596       REAL(wp)     ::  q_s               !<
2597       REAL(wp)     ::  sat               !<
2598       REAL(wp)     ::  t_l               !<
2599       REAL(wp)     ::  temp              !<
2600       REAL(wp)     ::  xr                !<
2601
2602       DO  k = nzb+1, nzt
2603!
2604!--       Predetermine flag to mask topography
2605          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2606
2607          IF ( qr_1d(k) > eps_sb )  THEN
2608!
2609!--          Actual liquid water temperature:
2610             t_l = t_d_pt(k) * pt_1d(k)
2611!
2612!--          Saturation vapor pressure at t_l:
2613             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
2614                                    ( t_l - 35.86_wp )                         &
2615                                  )
2616!
2617!--          Computation of saturation humidity:
2618             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2619             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2620             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
2621!
2622!--          Supersaturation:
2623             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2624!
2625!--          Evaporation needs only to be calculated in subsaturated regions
2626             IF ( sat < 0.0_wp )  THEN
2627!
2628!--             Actual temperature:
2629                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
2630
2631                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
2632                                    ( thermal_conductivity_l * temp ) +        &
2633                                    r_v * temp / ( diff_coeff_l * e_s )        &
2634                                  )
2635!
2636!--             Mean weight of rain drops
2637                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
2638!
2639!--             Weight averaged diameter of rain drops:
2640                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2641!
2642!--             Compute ventilation factor and intercept parameter
2643!--             (Seifert and Beheng, 2006; Seifert, 2008):
2644                IF ( ventilation_effect )  THEN
2645!
2646!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2647!--                Stevens and Seifert, 2008):
2648                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2649!
2650!--                Slope parameter of gamma distribution (Seifert, 2008):
2651                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
2652                                ( mu_r + 1.0_wp )                              &
2653                              )**( 1.0_wp / 3.0_wp ) / dr
2654
2655                   mu_r_2   = mu_r + 2.0_wp
2656                   mu_r_5d2 = mu_r + 2.5_wp
2657
2658                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
2659                            b_vent * schmidt_p_1d3 *                           &
2660                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
2661                            lambda_r**( -mu_r_5d2 ) *                          &
2662                            ( 1.0_wp -                                         &
2663                              0.5_wp * ( b_term / a_term ) *                   &
2664                              ( lambda_r / ( c_term + lambda_r )               &
2665                              )**mu_r_5d2 -                                    &
2666                              0.125_wp * ( b_term / a_term )**2 *              &
2667                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
2668                              )**mu_r_5d2 -                                    &
2669                              0.0625_wp * ( b_term / a_term )**3 *             &
2670                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
2671                              )**mu_r_5d2 -                                    &
2672                              0.0390625_wp * ( b_term / a_term )**4 *          &
2673                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
2674                              )**mu_r_5d2                                      &
2675                            )
2676
2677                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
2678                            gamm( mu_r + 1.0_wp )
2679                ELSE
2680                   f_vent = 1.0_wp
2681                   nr_0   = nr_1d(k) * dr
2682                ENDIF
2683!
2684!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
2685                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
2686                evap    = MAX( evap, -qr_1d(k) / dt_micro )
2687                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
2688                               -nr_1d(k) / dt_micro )
2689
2690                qr_1d(k) = qr_1d(k) + evap    * dt_micro * flag
2691                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
2692
2693             ENDIF
2694          ENDIF
2695
2696       ENDDO
2697
2698    END SUBROUTINE evaporation_rain_ij
2699
2700
2701!------------------------------------------------------------------------------!
2702! Description:
2703! ------------
2704!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2705!> Call for grid point i,j
2706!------------------------------------------------------------------------------!
2707    SUBROUTINE sedimentation_cloud_ij( i, j )
2708
2709       USE arrays_3d,                                                          &
2710           ONLY:  ddzu, dzu, prr
2711
2712       USE cloud_parameters,                                                   &
2713           ONLY:  hyrho, l_d_cp, pt_d_t
2714
2715       USE control_parameters,                                                 &
2716           ONLY:  call_microphysics_at_all_substeps,                           &
2717                  intermediate_timestep_count, microphysics_morrison
2718
2719       USE indices,                                                            &
2720           ONLY:  nzb, nzb, nzt, wall_flags_0
2721
2722       USE kinds
2723
2724       USE statistics,                                                         &
2725           ONLY:  weight_substep
2726
2727       IMPLICIT NONE
2728
2729       INTEGER(iwp) ::  i             !<
2730       INTEGER(iwp) ::  j             !<
2731       INTEGER(iwp) ::  k             !<
2732
2733       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
2734       REAL(wp)     ::  nc_sedi !<
2735
2736       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
2737       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
2738
2739       sed_qc(nzt+1) = 0.0_wp
2740       sed_nc(nzt+1) = 0.0_wp
2741
2742
2743       DO  k = nzt, nzb+1, -1
2744!
2745!--       Predetermine flag to mask topography
2746          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2747          nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
2748!
2749!--       Sedimentation fluxes for number concentration are only calculated
2750!--       for cloud_scheme = 'morrison'
2751          IF ( microphysics_morrison ) THEN
2752             IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
2753                sed_nc(k) = sed_qc_const *                                     &
2754                            ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
2755                            ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
2756             ELSE
2757                sed_nc(k) = 0.0_wp
2758             ENDIF
2759
2760             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
2761                              nc_1d(k) / dt_micro + sed_nc(k+1)                &
2762                            ) * flag
2763
2764             nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) *               &
2765                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
2766          ENDIF
2767
2768          IF ( qc_1d(k) > eps_sb  .AND.  nc_sedi > eps_mr )  THEN
2769             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
2770                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
2771          ELSE
2772             sed_qc(k) = 0.0_wp
2773          ENDIF
2774
2775          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
2776                                      dt_micro + sed_qc(k+1)                   &
2777                         ) * flag
2778
2779          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2780                                hyrho(k) * dt_micro * flag
2781          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2782                                hyrho(k) * dt_micro * flag
2783          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2784                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
2785
2786!
2787!--       Compute the precipitation rate of cloud (fog) droplets
2788          IF ( call_microphysics_at_all_substeps )  THEN
2789             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
2790                              weight_substep(intermediate_timestep_count) * flag
2791          ELSE
2792             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2793          ENDIF
2794
2795       ENDDO
2796
2797    END SUBROUTINE sedimentation_cloud_ij
2798
2799
2800!------------------------------------------------------------------------------!
2801! Description:
2802! ------------
2803!> Computation of sedimentation flux. Implementation according to Stevens
2804!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
2805!------------------------------------------------------------------------------!
2806    SUBROUTINE sedimentation_rain_ij( i, j )
2807
2808       USE arrays_3d,                                                          &
2809           ONLY:  ddzu, dzu, prr
2810
2811       USE cloud_parameters,                                                   &
2812           ONLY:  hyrho, l_d_cp, pt_d_t
2813
2814       USE control_parameters,                                                 &
2815           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
2816
2817       USE indices,                                                            &
2818           ONLY:  nzb, nzb, nzt, wall_flags_0
2819
2820       USE kinds
2821
2822       USE statistics,                                                         &
2823           ONLY:  weight_substep
2824
2825       USE surface_mod,                                                        &
2826           ONLY :  bc_h
2827
2828       IMPLICIT NONE
2829
2830       INTEGER(iwp) ::  i             !< running index x direction
2831       INTEGER(iwp) ::  j             !< running index y direction
2832       INTEGER(iwp) ::  k             !< running index z direction
2833       INTEGER(iwp) ::  k_run         !<
2834       INTEGER(iwp) ::  m             !< running index surface elements
2835       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2836       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2837
2838       REAL(wp)     ::  c_run                      !<
2839       REAL(wp)     ::  d_max                      !<
2840       REAL(wp)     ::  d_mean                     !<
2841       REAL(wp)     ::  d_min                      !<
2842       REAL(wp)     ::  dr                         !<
2843       REAL(wp)     ::  flux                       !<
2844       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
2845       REAL(wp)     ::  lambda_r                   !<
2846       REAL(wp)     ::  mu_r                       !<
2847       REAL(wp)     ::  z_run                      !<
2848
2849       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2850       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2851       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2852       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2853       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2854       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2855       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2856       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
2857
2858!
2859!--    Compute velocities
2860       DO  k = nzb+1, nzt
2861!
2862!--       Predetermine flag to mask topography
2863          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2864
2865          IF ( qr_1d(k) > eps_sb )  THEN
2866!
2867!--          Weight averaged diameter of rain drops:
2868             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
2869!
2870!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2871!--          Stevens and Seifert, 2008):
2872             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2873!
2874!--          Slope parameter of gamma distribution (Seifert, 2008):
2875             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
2876                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2877
2878             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2879                                         a_term - b_term * ( 1.0_wp +          &
2880                                            c_term / lambda_r )**( -1.0_wp *   &
2881                                            ( mu_r + 1.0_wp ) )                &
2882                                        )                                      &
2883                          ) * flag
2884             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2885                                         a_term - b_term * ( 1.0_wp +          &
2886                                            c_term / lambda_r )**( -1.0_wp *   &
2887                                            ( mu_r + 4.0_wp ) )                &
2888                                       )                                       &
2889                          ) * flag
2890          ELSE
2891             w_nr(k) = 0.0_wp
2892             w_qr(k) = 0.0_wp
2893          ENDIF
2894       ENDDO
2895!
2896!--    Adjust boundary values using surface data type.
2897!--    Upward facing non-natural
2898       surf_s = bc_h(0)%start_index(j,i)
2899       surf_e = bc_h(0)%end_index(j,i)
2900       DO  m = surf_s, surf_e
2901          k         = bc_h(0)%k(m)
2902          w_nr(k-1) = w_nr(k)
2903          w_qr(k-1) = w_qr(k)
2904       ENDDO
2905!
2906!--    Downward facing non-natural
2907       surf_s = bc_h(1)%start_index(j,i)
2908       surf_e = bc_h(1)%end_index(j,i)
2909       DO  m = surf_s, surf_e
2910          k         = bc_h(1)%k(m)
2911          w_nr(k+1) = w_nr(k)
2912          w_qr(k+1) = w_qr(k)
2913       ENDDO
2914!
2915!--    Neumann boundary condition at model top
2916       w_nr(nzt+1) = 0.0_wp
2917       w_qr(nzt+1) = 0.0_wp
2918!
2919!--    Compute Courant number
2920       DO  k = nzb+1, nzt
2921!
2922!--       Predetermine flag to mask topography
2923          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2924
2925          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
2926                    dt_micro * ddzu(k) * flag
2927          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
2928                    dt_micro * ddzu(k) * flag
2929       ENDDO
2930!
2931!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2932       IF ( limiter_sedimentation )  THEN
2933
2934          DO k = nzb+1, nzt
2935!
2936!--          Predetermine flag to mask topography
2937             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2938
2939             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
2940             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2941             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
2942
2943             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2944                                                        2.0_wp * d_max,        &
2945                                                        ABS( d_mean ) ) * flag
2946
2947             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
2948             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2949             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
2950
2951             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2952                                                        2.0_wp * d_max,        &
2953                                                        ABS( d_mean ) ) * flag
2954          ENDDO
2955
2956       ELSE
2957
2958          nr_slope = 0.0_wp
2959          qr_slope = 0.0_wp
2960
2961       ENDIF
2962
2963       sed_nr(nzt+1) = 0.0_wp
2964       sed_qr(nzt+1) = 0.0_wp
2965!
2966!--    Compute sedimentation flux
2967       DO  k = nzt, nzb+1, -1
2968!
2969!--       Predetermine flag to mask topography
2970          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2971!
2972!--       Sum up all rain drop number densities which contribute to the flux
2973!--       through k-1/2
2974          flux  = 0.0_wp
2975          z_run = 0.0_wp ! height above z(k)
2976          k_run = k
2977          c_run = MIN( 1.0_wp, c_nr(k) )
2978          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
2979             flux  = flux + hyrho(k_run) *                                     &
2980                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
2981                     0.5_wp ) * c_run * dzu(k_run) * flag
2982             z_run = z_run + dzu(k_run)            * flag
2983             k_run = k_run + 1                     * flag
2984             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
2985          ENDDO
2986!
2987!--       It is not allowed to sediment more rain drop number density than
2988!--       available
2989          flux = MIN( flux,                                                    &
2990                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
2991
2992          sed_nr(k) = flux / dt_micro * flag
2993          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
2994                                    hyrho(k) * dt_micro * flag
2995!
2996!--       Sum up all rain water content which contributes to the flux
2997!--       through k-1/2
2998          flux  = 0.0_wp
2999          z_run = 0.0_wp ! height above z(k)
3000          k_run = k
3001          c_run = MIN( 1.0_wp, c_qr(k) )
3002
3003          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3004
3005             flux  = flux + hyrho(k_run) *                                     &
3006                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
3007                     0.5_wp ) * c_run * dzu(k_run) * flag
3008             z_run = z_run + dzu(k_run)            * flag
3009             k_run = k_run + 1                     * flag
3010             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
3011
3012          ENDDO
3013!
3014!--       It is not allowed to sediment more rain water content than available
3015          flux = MIN( flux,                                                    &
3016                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
3017
3018          sed_qr(k) = flux / dt_micro * flag
3019
3020          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
3021                                hyrho(k) * dt_micro * flag
3022          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
3023                                hyrho(k) * dt_micro * flag
3024          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
3025                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
3026!
3027!--       Compute the rain rate
3028          IF ( call_microphysics_at_all_substeps )  THEN
3029             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
3030                          * weight_substep(intermediate_timestep_count) * flag
3031          ELSE
3032             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3033          ENDIF
3034
3035       ENDDO
3036
3037    END SUBROUTINE sedimentation_rain_ij
3038
3039
3040!------------------------------------------------------------------------------!
3041! Description:
3042! ------------
3043!> This subroutine computes the precipitation amount due to gravitational
3044!> settling of rain and cloud (fog) droplets
3045!------------------------------------------------------------------------------!
3046    SUBROUTINE calc_precipitation_amount_ij( i, j )
3047
3048       USE arrays_3d,                                                          &
3049           ONLY:  precipitation_amount, prr
3050
3051       USE cloud_parameters,                                                   &
3052           ONLY:  hyrho
3053
3054       USE control_parameters,                                                 &
3055           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
3056                  intermediate_timestep_count, intermediate_timestep_count_max,&
3057                  precipitation_amount_interval, time_do2d_xy
3058
3059       USE indices,                                                            &
3060           ONLY:  nzb, nzt, wall_flags_0
3061
3062       USE kinds
3063
3064       USE surface_mod,                                                        &
3065           ONLY :  bc_h
3066
3067       IMPLICIT NONE
3068
3069       INTEGER(iwp) ::  i             !< running index x direction
3070       INTEGER(iwp) ::  j             !< running index y direction
3071       INTEGER(iwp) ::  k             !< running index z direction
3072       INTEGER(iwp) ::  m             !< running index surface elements
3073       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
3074       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
3075
3076       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3077            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
3078            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3079       THEN
3080
3081          surf_s = bc_h(0)%start_index(j,i)
3082          surf_e = bc_h(0)%end_index(j,i)
3083          DO  m = surf_s, surf_e
3084             k                         = bc_h(0)%k(m)
3085             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
3086                                               prr(k,j,i) * hyrho(k) * dt_3d
3087          ENDDO
3088
3089       ENDIF
3090
3091    END SUBROUTINE calc_precipitation_amount_ij
3092
3093!------------------------------------------------------------------------------!
3094! Description:
3095! ------------
3096!> This function computes the gamma function (Press et al., 1992).
3097!> The gamma function is needed for the calculation of the evaporation
3098!> of rain drops.
3099!------------------------------------------------------------------------------!
3100    FUNCTION gamm( xx )
3101
3102       USE kinds
3103
3104       IMPLICIT NONE
3105
3106       INTEGER(iwp) ::  j            !<
3107
3108       REAL(wp)     ::  gamm         !<
3109       REAL(wp)     ::  ser          !<
3110       REAL(wp)     ::  tmp          !<
3111       REAL(wp)     ::  x_gamm       !<
3112       REAL(wp)     ::  xx           !<
3113       REAL(wp)     ::  y_gamm       !<
3114
3115
3116       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
3117       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
3118                                           -86.50532032941677_wp,      &
3119                                            24.01409824083091_wp,      &
3120                                            -1.231739572450155_wp,     &
3121                                             0.1208650973866179E-2_wp, &
3122                                            -0.5395239384953E-5_wp /)     !<
3123
3124       x_gamm = xx
3125       y_gamm = x_gamm
3126       tmp = x_gamm + 5.5_wp
3127       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
3128       ser = 1.000000000190015_wp
3129
3130       DO  j = 1, 6
3131          y_gamm = y_gamm + 1.0_wp
3132          ser    = ser + cof( j ) / y_gamm
3133       ENDDO
3134
3135!
3136!--    Until this point the algorithm computes the logarithm of the gamma
3137!--    function. Hence, the exponential function is used.
3138!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
3139       gamm = EXP( tmp ) * stp * ser / x_gamm
3140
3141       RETURN
3142
3143    END FUNCTION gamm
3144
3145 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.