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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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