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

Last change on this file since 3163 was 3026, checked in by schwenkel, 6 years ago

Changed the name specific humidity to mixing ratio

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