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

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

Merge of branch palm4u into trunk

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