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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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