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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

  • Property svn:keywords set to Id
File size: 95.9 KB
Line 
1!> @file microphysics_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Adjustments to new topography and surface concept
23!
24! Former revisions:
25! -----------------
26! $Id: microphysics_mod.f90 2232 2017-05-30 17:47:52Z suehring $
27!
28! 2155 2017-02-21 09:57:40Z hoffmann
29! Bugfix in the calculation of microphysical quantities on ghost points.
30!
31! 2031 2016-10-21 15:11:58Z knoop
32! renamed variable rho to rho_ocean
33!
34! 2000 2016-08-20 18:09:15Z knoop
35! Forced header and separation lines into 80 columns
36!
37! 1850 2016-04-08 13:29:27Z maronga
38! Module renamed
39! Adapted for modularization of microphysics.
40!
41! 1845 2016-04-08 08:29:13Z raasch
42! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
43! point (instead of one point above surface)
44!
45! 1831 2016-04-07 13:15:51Z hoffmann
46! turbulence renamed collision_turbulence,
47! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
48! avaialble for microphysics_kessler.
49!
50! 1822 2016-04-07 07:49:42Z hoffmann
51! Unused variables removed.
52! Kessler scheme integrated.
53!
54! 1691 2015-10-26 16:17:44Z maronga
55! Added new routine calc_precipitation_amount. The routine now allows to account
56! for precipitation due to sedimenation of cloud (fog) droplets
57!
58! 1682 2015-10-07 23:56:08Z knoop
59! Code annotations made doxygen readable
60!
61! 1646 2015-09-02 16:00:10Z hoffmann
62! Bugfix: Wrong computation of d_mean.
63!
64! 1361 2014-04-16 15:17:48Z hoffmann
65! Bugfix in sedimentation_rain: Index corrected.
66! Vectorized version of adjust_cloud added.
67! Little reformatting of the code.
68!
69! 1353 2014-04-08 15:21:23Z heinze
70! REAL constants provided with KIND-attribute
71!
72! 1346 2014-03-27 13:18:20Z heinze
73! Bugfix: REAL constants provided with KIND-attribute especially in call of
74! intrinsic function like MAX, MIN, SIGN
75!
76! 1334 2014-03-25 12:21:40Z heinze
77! Bugfix: REAL constants provided with KIND-attribute
78!
79! 1322 2014-03-20 16:38:49Z raasch
80! REAL constants defined as wp-kind
81!
82! 1320 2014-03-20 08:40:49Z raasch
83! ONLY-attribute added to USE-statements,
84! kind-parameters added to all INTEGER and REAL declaration statements,
85! kinds are defined in new module kinds,
86! comment fields (!:) to be used for variable explanations added to
87! all variable declaration statements
88!
89! 1241 2013-10-30 11:36:58Z heinze
90! hyp and rho_ocean have to be calculated at each time step if data from external
91! file LSF_DATA are used
92!
93! 1115 2013-03-26 18:16:16Z hoffmann
94! microphyical tendencies are calculated in microphysics_control in an optimized
95! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
96!
97! 1106 2013-03-04 05:31:38Z raasch
98! small changes in code formatting
99!
100! 1092 2013-02-02 11:24:22Z raasch
101! unused variables removed
102! file put under GPL
103!
104! 1065 2012-11-22 17:42:36Z hoffmann
105! Sedimentation process implemented according to Stevens and Seifert (2008).
106! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
107! and Stevens, 2010).
108!
109! 1053 2012-11-13 17:11:03Z hoffmann
110! initial revision
111!
112! Description:
113! ------------
114!> Calculate bilk cloud microphysics.
115!------------------------------------------------------------------------------!
116 MODULE microphysics_mod
117
118    USE  kinds
119
120    IMPLICIT NONE
121
122    LOGICAL ::  cloud_water_sedimentation = .FALSE.  !< cloud water sedimentation
123    LOGICAL ::  limiter_sedimentation = .TRUE.       !< sedimentation limiter
124    LOGICAL ::  collision_turbulence = .FALSE.       !< turbulence effects
125    LOGICAL ::  ventilation_effect = .TRUE.          !< ventilation effect
126
127    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
128    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
129    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
130    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
131    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
132    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
133    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
134    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
135    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
136    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
137    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
138    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
139    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
140    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
141    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
142    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
143    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
144    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
145    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
146    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
147    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
148    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
149    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
150    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
151    REAL(wp) ::  kappa_rr = 60.7_wp        !< const. in collision kernel (kg-1/3)
152    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
153    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
154    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
155    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
156    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
157    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
158    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
159    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
160    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
161    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
162
163    REAL(wp) ::  c_sedimentation = 2.0_wp  !< Courant number of sedimentation process
164    REAL(wp) ::  dpirho_l                  !< 6.0 / ( pi * rho_l )
165    REAL(wp) ::  dt_micro                  !< microphysics time step
166    REAL(wp) ::  nc_const = 70.0E6_wp      !< cloud droplet concentration
167    REAL(wp) ::  dt_precipitation = 100.0_wp !< timestep precipitation (s)
168    REAL(wp) ::  sed_qc_const              !< const. for sedimentation of cloud water
169    REAL(wp) ::  pirho_l                   !< pi * rho_l / 6.0;
170
171    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nc_1d  !<
172    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nr_1d  !<
173    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_1d  !<
174    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qc_1d  !<
175    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qr_1d  !<
176    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_1d   !<
177
178    SAVE
179
180    PRIVATE
181    PUBLIC microphysics_control, microphysics_init
182
183    PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation,   &
184           dt_precipitation, limiter_sedimentation, nc_const, sigma_gc,        &
185           ventilation_effect
186
187    INTERFACE microphysics_control
188       MODULE PROCEDURE microphysics_control
189       MODULE PROCEDURE microphysics_control_ij
190    END INTERFACE microphysics_control
191
192    INTERFACE adjust_cloud
193       MODULE PROCEDURE adjust_cloud
194       MODULE PROCEDURE adjust_cloud_ij
195    END INTERFACE adjust_cloud
196
197    INTERFACE autoconversion
198       MODULE PROCEDURE autoconversion
199       MODULE PROCEDURE autoconversion_ij
200    END INTERFACE autoconversion
201
202    INTERFACE autoconversion_kessler
203       MODULE PROCEDURE autoconversion_kessler
204       MODULE PROCEDURE autoconversion_kessler_ij
205    END INTERFACE autoconversion_kessler
206
207    INTERFACE accretion
208       MODULE PROCEDURE accretion
209       MODULE PROCEDURE accretion_ij
210    END INTERFACE accretion
211
212    INTERFACE selfcollection_breakup
213       MODULE PROCEDURE selfcollection_breakup
214       MODULE PROCEDURE selfcollection_breakup_ij
215    END INTERFACE selfcollection_breakup
216
217    INTERFACE evaporation_rain
218       MODULE PROCEDURE evaporation_rain
219       MODULE PROCEDURE evaporation_rain_ij
220    END INTERFACE evaporation_rain
221
222    INTERFACE sedimentation_cloud
223       MODULE PROCEDURE sedimentation_cloud
224       MODULE PROCEDURE sedimentation_cloud_ij
225    END INTERFACE sedimentation_cloud
226
227    INTERFACE sedimentation_rain
228       MODULE PROCEDURE sedimentation_rain
229       MODULE PROCEDURE sedimentation_rain_ij
230    END INTERFACE sedimentation_rain
231
232    INTERFACE calc_precipitation_amount
233       MODULE PROCEDURE calc_precipitation_amount
234       MODULE PROCEDURE calc_precipitation_amount_ij
235    END INTERFACE calc_precipitation_amount
236
237 CONTAINS
238!------------------------------------------------------------------------------!
239! Description:
240! ------------
241!> Initialization of bulk microphysics
242!------------------------------------------------------------------------------!
243    SUBROUTINE microphysics_init
244
245       USE arrays_3d,                                                          &
246           ONLY:  dzu
247
248       USE constants,                                                          &
249           ONLY:  pi
250
251       USE cloud_parameters,                                                   &
252           ONLY:  rho_l
253
254       USE control_parameters,                                                 &
255           ONLY:  microphysics_seifert
256
257       USE indices,                                                            &
258           ONLY:  nzb, nzt
259
260       IMPLICIT NONE
261
262!
263!--    constant for the sedimentation of cloud water (2-moment cloud physics)
264       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
265                          )**( 2.0_wp / 3.0_wp ) *                             &
266                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
267
268!
269!--    Calculate timestep according to precipitation
270       IF ( microphysics_seifert )  THEN
271          dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
272                             w_precipitation
273       ENDIF
274
275!
276!--    Pre-calculate frequently calculated fractions of pi and rho_l
277       pirho_l  = pi * rho_l / 6.0_wp
278       dpirho_l = 1.0_wp / pirho_l
279
280!
281!--    Allocate 1D microphysics arrays
282       ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
283                  qc_1d(nzb:nzt+1) )
284
285       IF ( microphysics_seifert )  THEN
286          ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
287       ENDIF
288
289!
290!--    Initialze nc_1d with nc_const
291       nc_1d = nc_const
292
293    END SUBROUTINE microphysics_init
294
295
296!------------------------------------------------------------------------------!
297! Description:
298! ------------
299!> Control of microphysics for all grid points
300!------------------------------------------------------------------------------!
301    SUBROUTINE microphysics_control
302
303       USE arrays_3d,                                                          &
304           ONLY:  hyp, pt_init, prr, zu
305
306       USE cloud_parameters,                                                   &
307           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
308
309       USE control_parameters,                                                 &
310           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
311                  intermediate_timestep_count, large_scale_forcing,            &
312                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
313                  pt_surface, rho_surface,surface_pressure
314
315       USE indices,                                                            &
316           ONLY:  nzb, nzt
317
318       USE kinds
319
320       USE statistics,                                                         &
321           ONLY:  weight_pres
322
323       IMPLICIT NONE
324
325       INTEGER(iwp) ::  k                 !<
326
327       REAL(wp)     ::  t_surface         !<
328
329       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
330!
331!--       Calculate:
332!--       pt / t : ratio of potential and actual temperature (pt_d_t)
333!--       t / pt : ratio of actual and potential temperature (t_d_pt)
334!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
335          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
336          DO  k = nzb, nzt+1
337             hyp(k)    = surface_pressure * 100.0_wp * &
338                         ( ( t_surface - g / cp * zu(k) ) /                    &
339                         t_surface )**(1.0_wp / 0.286_wp)
340             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
341             t_d_pt(k) = 1.0_wp / pt_d_t(k)
342             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
343          ENDDO
344
345!
346!--       Compute reference density
347          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
348       ENDIF
349
350!
351!--    Compute length of time step
352       IF ( call_microphysics_at_all_substeps )  THEN
353          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
354       ELSE
355          dt_micro = dt_3d
356       ENDIF
357
358!
359!--    Reset precipitation rate
360       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
361
362!
363!--    Compute cloud physics
364       IF ( microphysics_kessler )  THEN
365
366          CALL autoconversion_kessler
367          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
368
369       ELSEIF ( microphysics_seifert )  THEN
370
371          CALL adjust_cloud
372          CALL autoconversion
373          CALL accretion
374          CALL selfcollection_breakup
375          CALL evaporation_rain
376          CALL sedimentation_rain
377          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
378
379       ENDIF
380
381       CALL calc_precipitation_amount
382
383    END SUBROUTINE microphysics_control
384
385!------------------------------------------------------------------------------!
386! Description:
387! ------------
388!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
389!> evaporation of rain drops due to too small or too big weights
390!> of rain drops (Stevens and Seifert, 2008).
391!------------------------------------------------------------------------------!
392    SUBROUTINE adjust_cloud
393
394       USE arrays_3d,                                                          &
395           ONLY:  qr, nr
396
397       USE cloud_parameters,                                                   &
398           ONLY:  hyrho
399
400       USE cpulog,                                                             &
401           ONLY:  cpu_log, log_point_s
402
403       USE indices,                                                            &
404           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
405
406       USE kinds
407
408       IMPLICIT NONE
409
410       INTEGER(iwp) ::  i                 !<
411       INTEGER(iwp) ::  j                 !<
412       INTEGER(iwp) ::  k                 !<
413
414       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
415
416       DO  i = nxlg, nxrg
417          DO  j = nysg, nyng
418             DO  k = nzb+1, nzt
419                IF ( qr(k,j,i) <= eps_sb )  THEN
420                   qr(k,j,i) = 0.0_wp
421                   nr(k,j,i) = 0.0_wp
422                ELSE
423                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
424                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
425                                       MERGE( 1.0_wp, 0.0_wp,                  &           
426                                              BTEST( wall_flags_0(k,j,i), 0 ) )
427                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
428                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
429                                       MERGE( 1.0_wp, 0.0_wp,                  &           
430                                              BTEST( wall_flags_0(k,j,i), 0 ) )
431                   ENDIF
432                ENDIF
433             ENDDO
434          ENDDO
435       ENDDO
436
437       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
438
439    END SUBROUTINE adjust_cloud
440
441
442!------------------------------------------------------------------------------!
443! Description:
444! ------------
445!> Autoconversion rate (Seifert and Beheng, 2006).
446!------------------------------------------------------------------------------!
447    SUBROUTINE autoconversion
448
449       USE arrays_3d,                                                          &
450           ONLY:  diss, dzu, nr, qc, qr
451
452       USE cloud_parameters,                                                   &
453           ONLY:  hyrho
454
455       USE control_parameters,                                                 &
456           ONLY:  rho_surface
457
458       USE cpulog,                                                             &
459           ONLY:  cpu_log, log_point_s
460
461       USE grid_variables,                                                     &
462           ONLY:  dx, dy
463
464       USE indices,                                                            &
465           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
466
467       USE kinds
468
469       IMPLICIT NONE
470
471       INTEGER(iwp) ::  i                 !<
472       INTEGER(iwp) ::  j                 !<
473       INTEGER(iwp) ::  k                 !<
474
475       REAL(wp)     ::  alpha_cc          !<
476       REAL(wp)     ::  autocon           !<
477       REAL(wp)     ::  dissipation       !<
478       REAL(wp)     ::  flag              !< flag to mask topography grid points
479       REAL(wp)     ::  k_au              !<
480       REAL(wp)     ::  l_mix             !<
481       REAL(wp)     ::  nu_c              !<
482       REAL(wp)     ::  phi_au            !<
483       REAL(wp)     ::  r_cc              !<
484       REAL(wp)     ::  rc                !<
485       REAL(wp)     ::  re_lambda         !<
486       REAL(wp)     ::  sigma_cc          !<
487       REAL(wp)     ::  tau_cloud         !<
488       REAL(wp)     ::  xc                !<
489
490       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
491
492       DO  i = nxlg, nxrg
493          DO  j = nysg, nyng
494             DO  k = nzb+1, nzt
495!
496!--             Predetermine flag to mask topography
497                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
498
499                IF ( qc(k,j,i) > eps_sb )  THEN
500
501                   k_au = k_cc / ( 20.0_wp * x0 )
502!
503!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
504!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
505                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
506!
507!--                Universal function for autoconversion process
508!--                (Seifert and Beheng, 2006):
509                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
510                            ( 1.0_wp - tau_cloud**0.68_wp )**3
511!
512!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
513!--                (Use constant nu_c = 1.0_wp instead?)
514                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
515!
516!--                Mean weight of cloud droplets:
517                   xc = hyrho(k) * qc(k,j,i) / nc_const
518!
519!--                Parameterized turbulence effects on autoconversion (Seifert,
520!--                Nuijens and Stevens, 2010)
521                   IF ( collision_turbulence )  THEN
522!
523!--                   Weight averaged radius of cloud droplets:
524                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
525
526                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
527                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
528                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
529!
530!--                   Mixing length (neglecting distance to ground and
531!--                   stratification)
532                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
533!
534!--                   Limit dissipation rate according to Seifert, Nuijens and
535!--                   Stevens (2010)
536                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
537!
538!--                   Compute Taylor-microscale Reynolds number:
539                      re_lambda = 6.0_wp / 11.0_wp *                           &
540                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
541                                  SQRT( 15.0_wp / kin_vis_air ) *              &
542                                  dissipation**( 1.0_wp / 6.0_wp )
543!
544!--                   The factor of 1.0E4 is needed to convert the dissipation
545!--                   rate from m2 s-3 to cm2 s-3.
546                      k_au = k_au * ( 1.0_wp +                                 &
547                                      dissipation * 1.0E4_wp *                 &
548                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
549                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
550                                                                      r_cc ) / &
551                                                        sigma_cc )**2          &
552                                                      ) + beta_cc              &
553                                      )                                        &
554                                    )
555                   ENDIF
556!
557!--                Autoconversion rate (Seifert and Beheng, 2006):
558                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
559                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
560                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
561                             rho_surface
562                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
563
564                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
565                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
566                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
567                                                              * flag
568
569                ENDIF
570
571             ENDDO
572          ENDDO
573       ENDDO
574
575       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
576
577    END SUBROUTINE autoconversion
578
579
580!------------------------------------------------------------------------------!
581! Description:
582! ------------
583!> Autoconversion process (Kessler, 1969).
584!------------------------------------------------------------------------------!
585    SUBROUTINE autoconversion_kessler
586
587       USE arrays_3d,                                                          &
588           ONLY:  dzw, pt, prr, q, qc
589
590       USE cloud_parameters,                                                   &
591           ONLY:  l_d_cp, pt_d_t
592
593       USE indices,                                                            &
594           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
595
596       USE kinds
597
598
599       IMPLICIT NONE
600
601       INTEGER(iwp) ::  i      !<
602       INTEGER(iwp) ::  j      !<
603       INTEGER(iwp) ::  k      !<
604       INTEGER(iwp) ::  k_wall !< topgraphy top index
605
606       REAL(wp)    ::  dqdt_precip !<
607       REAL(wp)    ::  flag        !< flag to mask topography grid points
608
609       DO  i = nxlg, nxrg
610          DO  j = nysg, nyng
611!
612!--          Determine vertical index of topography top
613             k_wall = MAXLOC(                                                  &
614                          MERGE( 1, 0,                                         &
615                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
616                               ), DIM = 1                                      &
617                            ) - 1
618             DO  k = nzb+1, nzt
619!
620!--             Predetermine flag to mask topography
621                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
622
623                IF ( qc(k,j,i) > ql_crit )  THEN
624                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
625                ELSE
626                   dqdt_precip = 0.0_wp
627                ENDIF
628
629                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
630                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
631                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
632                                        pt_d_t(k)              * flag
633
634!
635!--             Compute the rain rate (stored on surface grid point)
636                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
637
638             ENDDO
639          ENDDO
640       ENDDO
641
642   END SUBROUTINE autoconversion_kessler
643
644
645!------------------------------------------------------------------------------!
646! Description:
647! ------------
648!> Accretion rate (Seifert and Beheng, 2006).
649!------------------------------------------------------------------------------!
650    SUBROUTINE accretion
651
652       USE arrays_3d,                                                          &
653           ONLY:  diss, qc, qr
654
655       USE cloud_parameters,                                                   &
656           ONLY:  hyrho
657
658       USE control_parameters,                                                 &
659           ONLY:  rho_surface
660
661       USE cpulog,                                                             &
662           ONLY:  cpu_log, log_point_s
663
664       USE indices,                                                            &
665           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
666
667       USE kinds
668
669       IMPLICIT NONE
670
671       INTEGER(iwp) ::  i                 !<
672       INTEGER(iwp) ::  j                 !<
673       INTEGER(iwp) ::  k                 !<
674
675       REAL(wp)     ::  accr              !<
676       REAL(wp)     ::  flag              !< flag to mask topography grid points
677       REAL(wp)     ::  k_cr              !<
678       REAL(wp)     ::  phi_ac            !<
679       REAL(wp)     ::  tau_cloud         !<
680
681       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
682
683       DO  i = nxlg, nxrg
684          DO  j = nysg, nyng
685             DO  k = nzb+1, nzt
686!
687!--             Predetermine flag to mask topography
688                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
689
690                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
691!
692!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
693                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
694!
695!--                Universal function for accretion process (Seifert and
696!--                Beheng, 2001):
697                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
698!
699!--                Parameterized turbulence effects on autoconversion (Seifert,
700!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
701!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
702                   IF ( collision_turbulence )  THEN
703                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
704                                       MIN( 600.0_wp,                          &
705                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
706                                     )
707                   ELSE
708                      k_cr = k_cr0
709                   ENDIF
710!
711!--                Accretion rate (Seifert and Beheng, 2006):
712                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
713                          SQRT( rho_surface * hyrho(k) )
714                   accr = MIN( accr, qc(k,j,i) / dt_micro )
715
716                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
717                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
718
719                ENDIF
720
721             ENDDO
722          ENDDO
723       ENDDO
724
725       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
726
727    END SUBROUTINE accretion
728
729
730!------------------------------------------------------------------------------!
731! Description:
732! ------------
733!> Collisional breakup rate (Seifert, 2008).
734!------------------------------------------------------------------------------!
735    SUBROUTINE selfcollection_breakup
736
737       USE arrays_3d,                                                          &
738           ONLY:  nr, qr
739
740       USE cloud_parameters,                                                   &
741           ONLY:  hyrho
742
743       USE control_parameters,                                                 &
744           ONLY:  rho_surface
745
746       USE cpulog,                                                             &
747           ONLY:  cpu_log, log_point_s
748
749       USE indices,                                                            &
750           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
751
752       USE kinds
753
754       IMPLICIT NONE
755
756       INTEGER(iwp) ::  i                 !<
757       INTEGER(iwp) ::  j                 !<
758       INTEGER(iwp) ::  k                 !<
759
760       REAL(wp)     ::  breakup           !<
761       REAL(wp)     ::  dr                !<
762       REAL(wp)     ::  flag              !< flag to mask topography grid points
763       REAL(wp)     ::  phi_br            !<
764       REAL(wp)     ::  selfcoll          !<
765
766       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
767
768       DO  i = nxlg, nxrg
769          DO  j = nysg, nyng
770             DO  k = nzb+1, nzt
771!
772!--             Predetermine flag to mask topography
773                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
774
775                IF ( qr(k,j,i) > eps_sb )  THEN
776!
777!--                Selfcollection rate (Seifert and Beheng, 2001):
778                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
779                              SQRT( hyrho(k) * rho_surface )
780!
781!--                Weight averaged diameter of rain drops:
782                   dr = ( hyrho(k) * qr(k,j,i) /                               &
783                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
784!
785!--                Collisional breakup rate (Seifert, 2008):
786                   IF ( dr >= 0.3E-3_wp )  THEN
787                      phi_br  = k_br * ( dr - 1.1E-3_wp )
788                      breakup = selfcoll * ( phi_br + 1.0_wp )
789                   ELSE
790                      breakup = 0.0_wp
791                   ENDIF
792
793                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
794                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
795
796                ENDIF
797             ENDDO
798          ENDDO
799       ENDDO
800
801       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
802
803    END SUBROUTINE selfcollection_breakup
804
805
806!------------------------------------------------------------------------------!
807! Description:
808! ------------
809!> Evaporation of precipitable water. Condensation is neglected for
810!> precipitable water.
811!------------------------------------------------------------------------------!
812    SUBROUTINE evaporation_rain
813
814       USE arrays_3d,                                                          &
815           ONLY:  hyp, nr, pt, q,  qc, qr
816
817       USE cloud_parameters,                                                   &
818           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
819
820       USE constants,                                                          &
821           ONLY:  pi
822
823       USE cpulog,                                                             &
824           ONLY:  cpu_log, log_point_s
825
826       USE indices,                                                            &
827           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
828
829       USE kinds
830
831       IMPLICIT NONE
832
833       INTEGER(iwp) ::  i                 !<
834       INTEGER(iwp) ::  j                 !<
835       INTEGER(iwp) ::  k                 !<
836
837       REAL(wp)     ::  alpha             !<
838       REAL(wp)     ::  dr                !<
839       REAL(wp)     ::  e_s               !<
840       REAL(wp)     ::  evap              !<
841       REAL(wp)     ::  evap_nr           !<
842       REAL(wp)     ::  f_vent            !<
843       REAL(wp)     ::  flag              !< flag to mask topography grid points
844       REAL(wp)     ::  g_evap            !<
845       REAL(wp)     ::  lambda_r          !<
846       REAL(wp)     ::  mu_r              !<
847       REAL(wp)     ::  mu_r_2            !<
848       REAL(wp)     ::  mu_r_5d2          !<
849       REAL(wp)     ::  nr_0              !<
850       REAL(wp)     ::  q_s               !<
851       REAL(wp)     ::  sat               !<
852       REAL(wp)     ::  t_l               !<
853       REAL(wp)     ::  temp              !<
854       REAL(wp)     ::  xr                !<
855
856       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
857
858       DO  i = nxlg, nxrg
859          DO  j = nysg, nyng
860             DO  k = nzb+1, nzt
861!
862!--             Predetermine flag to mask topography
863                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
864
865                IF ( qr(k,j,i) > eps_sb )  THEN
866!
867!--                Actual liquid water temperature:
868                   t_l = t_d_pt(k) * pt(k,j,i)
869!
870!--                Saturation vapor pressure at t_l:
871                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
872                                          ( t_l - 35.86_wp )                   &
873                                        )
874!
875!--                Computation of saturation humidity:
876                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
877                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
878                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
879                           ( 1.0_wp + alpha * q_s )
880!
881!--                Supersaturation:
882                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
883!
884!--                Evaporation needs only to be calculated in subsaturated regions
885                   IF ( sat < 0.0_wp )  THEN
886!
887!--                   Actual temperature:
888                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
889
890                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
891                                          l_v / ( thermal_conductivity_l * temp ) &
892                                          + r_v * temp / ( diff_coeff_l * e_s )   &
893                                        )
894!
895!--                   Mean weight of rain drops
896                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
897!
898!--                   Weight averaged diameter of rain drops:
899                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
900!
901!--                   Compute ventilation factor and intercept parameter
902!--                   (Seifert and Beheng, 2006; Seifert, 2008):
903                      IF ( ventilation_effect )  THEN
904!
905!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
906!--                      2005; Stevens and Seifert, 2008):
907                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
908                                                          ( dr - 1.4E-3_wp ) ) )
909!
910!--                      Slope parameter of gamma distribution (Seifert, 2008):
911                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
912                                      ( mu_r + 1.0_wp )                        &
913                                    )**( 1.0_wp / 3.0_wp ) / dr
914
915                         mu_r_2   = mu_r + 2.0_wp
916                         mu_r_5d2 = mu_r + 2.5_wp
917
918                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
919                                  lambda_r**( -mu_r_2 ) + b_vent *              &
920                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
921                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
922                                  ( 1.0_wp -                                    &
923                                    0.5_wp * ( b_term / a_term ) *              &
924                                    ( lambda_r / ( c_term + lambda_r )          &
925                                    )**mu_r_5d2 -                               &
926                                    0.125_wp * ( b_term / a_term )**2 *         &
927                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
928                                    )**mu_r_5d2 -                               &
929                                    0.0625_wp * ( b_term / a_term )**3 *        &
930                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
931                                    )**mu_r_5d2 -                               &
932                                    0.0390625_wp * ( b_term / a_term )**4 *     &
933                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
934                                    )**mu_r_5d2                                 &
935                                  )
936
937                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
938                                  gamm( mu_r + 1.0_wp )
939                      ELSE
940                         f_vent = 1.0_wp
941                         nr_0   = nr(k,j,i) * dr
942                      ENDIF
943   !
944   !--                Evaporation rate of rain water content (Seifert and
945   !--                Beheng, 2006):
946                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
947                                hyrho(k)
948                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
949                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
950                                     -nr(k,j,i) / dt_micro )
951
952                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
953                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
954
955                   ENDIF
956                ENDIF
957
958             ENDDO
959          ENDDO
960       ENDDO
961
962       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
963
964    END SUBROUTINE evaporation_rain
965
966
967!------------------------------------------------------------------------------!
968! Description:
969! ------------
970!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
971!------------------------------------------------------------------------------!
972    SUBROUTINE sedimentation_cloud
973
974       USE arrays_3d,                                                          &
975           ONLY:  ddzu, dzu, pt, prr, q, qc
976
977       USE cloud_parameters,                                                   &
978           ONLY:  hyrho, l_d_cp, pt_d_t
979
980       USE control_parameters,                                                 &
981           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
982
983       USE cpulog,                                                             &
984           ONLY:  cpu_log, log_point_s
985
986       USE indices,                                                            &
987           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
988
989       USE kinds
990
991       USE statistics,                                                         &
992           ONLY:  weight_substep
993
994
995       IMPLICIT NONE
996
997       INTEGER(iwp) ::  i             !<
998       INTEGER(iwp) ::  j             !<
999       INTEGER(iwp) ::  k             !<
1000
1001       REAL(wp)                       ::  flag   !< flag to mask topography grid points
1002       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
1003
1004       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
1005
1006       sed_qc(nzt+1) = 0.0_wp
1007
1008       DO  i = nxlg, nxrg
1009          DO  j = nysg, nyng
1010             DO  k = nzt, nzb+1, -1
1011!
1012!--             Predetermine flag to mask topography
1013                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1014
1015                IF ( qc(k,j,i) > eps_sb )  THEN
1016                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
1017                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
1018                                                                           flag
1019                ELSE
1020                   sed_qc(k) = 0.0_wp
1021                ENDIF
1022
1023                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
1024                                            dt_micro + sed_qc(k+1)             &
1025                               ) * flag
1026
1027                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
1028                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1029                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
1030                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1031                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
1032                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1033                                        pt_d_t(k) * dt_micro            * flag
1034
1035!
1036!--             Compute the precipitation rate due to cloud (fog) droplets
1037                IF ( call_microphysics_at_all_substeps )  THEN
1038                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
1039                                * weight_substep(intermediate_timestep_count)  &
1040                                * flag
1041                ELSE
1042                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
1043                ENDIF
1044
1045             ENDDO
1046          ENDDO
1047       ENDDO
1048
1049       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
1050
1051    END SUBROUTINE sedimentation_cloud
1052
1053
1054!------------------------------------------------------------------------------!
1055! Description:
1056! ------------
1057!> Computation of sedimentation flux. Implementation according to Stevens
1058!> and Seifert (2008). Code is based on UCLA-LES.
1059!------------------------------------------------------------------------------!
1060    SUBROUTINE sedimentation_rain
1061
1062       USE arrays_3d,                                                          &
1063           ONLY:  ddzu, dzu, nr, pt, prr, q, qr
1064
1065       USE cloud_parameters,                                                   &
1066           ONLY:  hyrho, l_d_cp, pt_d_t
1067
1068       USE control_parameters,                                                 &
1069           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
1070       USE cpulog,                                                             &
1071           ONLY:  cpu_log, log_point_s
1072
1073       USE indices,                                                            &
1074           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
1075
1076       USE kinds
1077
1078       USE statistics,                                                         &
1079           ONLY:  weight_substep
1080
1081       USE surface_mod,                                                        &
1082           ONLY :  bc_h
1083       
1084       IMPLICIT NONE
1085
1086       INTEGER(iwp) ::  i             !< running index x direction
1087       INTEGER(iwp) ::  j             !< running index y direction
1088       INTEGER(iwp) ::  k             !< running index z direction
1089       INTEGER(iwp) ::  k_run         !<
1090       INTEGER(iwp) ::  l             !< running index of surface type
1091       INTEGER(iwp) ::  m             !< running index surface elements
1092       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1093       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
1094
1095       REAL(wp)     ::  c_run                      !<
1096       REAL(wp)     ::  d_max                      !<
1097       REAL(wp)     ::  d_mean                     !<
1098       REAL(wp)     ::  d_min                      !<
1099       REAL(wp)     ::  dr                         !<
1100       REAL(wp)     ::  flux                       !<
1101       REAL(wp)     ::  flag                       !< flag to mask topography grid points
1102       REAL(wp)     ::  lambda_r                   !<
1103       REAL(wp)     ::  mu_r                       !<
1104       REAL(wp)     ::  z_run                      !<
1105
1106       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1107       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1108       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1109       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1110       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1111       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1112       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1113       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
1114
1115       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
1116
1117!
1118!--    Compute velocities
1119       DO  i = nxlg, nxrg
1120          DO  j = nysg, nyng
1121             DO  k = nzb+1, nzt
1122!
1123!--             Predetermine flag to mask topography
1124                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1125
1126                IF ( qr(k,j,i) > eps_sb )  THEN
1127!
1128!--                Weight averaged diameter of rain drops:
1129                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1130                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1131!
1132!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1133!--                Stevens and Seifert, 2008):
1134                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
1135                                                     ( dr - 1.4E-3_wp ) ) )
1136!
1137!--                Slope parameter of gamma distribution (Seifert, 2008):
1138                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1139                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1140
1141                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1142                                               a_term - b_term * ( 1.0_wp +    &
1143                                                  c_term /                     &
1144                                                  lambda_r )**( -1.0_wp *      &
1145                                                  ( mu_r + 1.0_wp ) )          &
1146                                              )                                &
1147                                ) * flag
1148
1149                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1150                                               a_term - b_term * ( 1.0_wp +    &
1151                                                  c_term /                     &
1152                                                  lambda_r )**( -1.0_wp *      &
1153                                                  ( mu_r + 4.0_wp ) )          &
1154                                             )                                 &
1155                                ) * flag
1156                ELSE
1157                   w_nr(k) = 0.0_wp
1158                   w_qr(k) = 0.0_wp
1159                ENDIF
1160             ENDDO
1161!
1162!--          Adjust boundary values using surface data type.
1163!--          Upward-facing
1164             surf_s = bc_h(0)%start_index(j,i)   
1165             surf_e = bc_h(0)%end_index(j,i)
1166             DO  m = surf_s, surf_e
1167                k         = bc_h(0)%k(m)
1168                w_nr(k-1) = w_nr(k)
1169                w_qr(k-1) = w_qr(k)
1170             ENDDO
1171!
1172!--          Downward-facing
1173             surf_s = bc_h(1)%start_index(j,i)   
1174             surf_e = bc_h(1)%end_index(j,i)
1175             DO  m = surf_s, surf_e
1176                k         = bc_h(1)%k(m)
1177                w_nr(k+1) = w_nr(k)
1178                w_qr(k+1) = w_qr(k)
1179             ENDDO
1180!
1181!--          Model top boundary value
1182             w_nr(nzt+1) = 0.0_wp
1183             w_qr(nzt+1) = 0.0_wp
1184!
1185!--          Compute Courant number
1186             DO  k = nzb+1, nzt
1187!
1188!--             Predetermine flag to mask topography
1189                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1190
1191                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
1192                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
1193                          dt_micro * ddzu(k) * flag
1194                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
1195                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
1196                          dt_micro * ddzu(k) * flag
1197             ENDDO
1198!
1199!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1200             IF ( limiter_sedimentation )  THEN
1201
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                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
1208                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1209                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1210
1211                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1212                                                              2.0_wp * d_max,  &
1213                                                              ABS( d_mean ) )  &
1214                                                      * flag
1215
1216                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
1217                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1218                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1219
1220                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1221                                                              2.0_wp * d_max,  &
1222                                                              ABS( d_mean ) )
1223                ENDDO
1224
1225             ELSE
1226
1227                nr_slope = 0.0_wp
1228                qr_slope = 0.0_wp
1229
1230             ENDIF
1231
1232             sed_nr(nzt+1) = 0.0_wp
1233             sed_qr(nzt+1) = 0.0_wp
1234!
1235!--          Compute sedimentation flux
1236             DO  k = nzt, nzb+1, -1
1237!
1238!--             Predetermine flag to mask topography
1239                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1240!
1241!--             Sum up all rain drop number densities which contribute to the flux
1242!--             through k-1/2
1243                flux  = 0.0_wp
1244                z_run = 0.0_wp ! height above z(k)
1245                k_run = k
1246                c_run = MIN( 1.0_wp, c_nr(k) )
1247                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1248                   flux  = flux + hyrho(k_run) *                               &
1249                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
1250                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
1251                                              * flag
1252                   z_run = z_run + dzu(k_run) * flag
1253                   k_run = k_run + 1          * flag
1254                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
1255                                              * flag
1256                ENDDO
1257!
1258!--             It is not allowed to sediment more rain drop number density than
1259!--             available
1260                flux = MIN( flux,                                              &
1261                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
1262                            dt_micro                                           &
1263                          )
1264
1265                sed_nr(k) = flux / dt_micro * flag
1266                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
1267                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1268!
1269!--             Sum up all rain water content which contributes to the flux
1270!--             through k-1/2
1271                flux  = 0.0_wp
1272                z_run = 0.0_wp ! height above z(k)
1273                k_run = k
1274                c_run = MIN( 1.0_wp, c_qr(k) )
1275
1276                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1277
1278                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
1279                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
1280                                  0.5_wp ) * c_run * dzu(k_run) * flag
1281                   z_run = z_run + dzu(k_run)                   * flag
1282                   k_run = k_run + 1                            * flag
1283                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
1284                                                                * flag
1285
1286                ENDDO
1287!
1288!--             It is not allowed to sediment more rain water content than
1289!--             available
1290                flux = MIN( flux,                                              &
1291                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
1292                            dt_micro                                           &
1293                          )
1294
1295                sed_qr(k) = flux / dt_micro * flag
1296
1297                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
1298                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1299                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
1300                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
1301                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
1302                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1303                                        pt_d_t(k) * dt_micro            * flag
1304!
1305!--             Compute the rain rate
1306                IF ( call_microphysics_at_all_substeps )  THEN
1307                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
1308                                * weight_substep(intermediate_timestep_count) &
1309                                * flag
1310                ELSE
1311                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
1312                ENDIF
1313
1314             ENDDO
1315          ENDDO
1316       ENDDO
1317
1318       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1319
1320    END SUBROUTINE sedimentation_rain
1321
1322
1323!------------------------------------------------------------------------------!
1324! Description:
1325! ------------
1326!> Computation of the precipitation amount due to gravitational settling of
1327!> rain and cloud (fog) droplets
1328!------------------------------------------------------------------------------!
1329    SUBROUTINE calc_precipitation_amount
1330
1331       USE arrays_3d,                                                          &
1332           ONLY:  precipitation_amount, prr
1333
1334       USE cloud_parameters,                                                   &
1335           ONLY:  hyrho
1336
1337       USE control_parameters,                                                 &
1338           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
1339                  intermediate_timestep_count, intermediate_timestep_count_max,&
1340                  precipitation_amount_interval, time_do2d_xy
1341
1342       USE indices,                                                            &
1343           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
1344
1345       USE kinds
1346
1347       USE surface_mod,                                                        &
1348           ONLY :  bc_h
1349
1350       IMPLICIT NONE
1351
1352       INTEGER(iwp) ::  i             !< running index x direction
1353       INTEGER(iwp) ::  j             !< running index y direction
1354       INTEGER(iwp) ::  k             !< running index y direction
1355       INTEGER(iwp) ::  m             !< running index surface elements
1356       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1357       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
1358
1359       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1360            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
1361            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1362       THEN
1363!
1364!--       Run over all upward-facing surface elements, i.e. non-natural,
1365!--       natural and urban
1366          DO  m = 1, bc_h(0)%ns
1367             i = bc_h(0)%i(m)           
1368             j = bc_h(0)%j(m)
1369             k = bc_h(0)%k(m)
1370             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
1371                                               prr(k,j,i) * hyrho(k) * dt_3d
1372          ENDDO
1373
1374       ENDIF
1375
1376    END SUBROUTINE calc_precipitation_amount
1377
1378
1379!------------------------------------------------------------------------------!
1380! Description:
1381! ------------
1382!> Control of microphysics for grid points i,j
1383!------------------------------------------------------------------------------!
1384
1385    SUBROUTINE microphysics_control_ij( i, j )
1386
1387       USE arrays_3d,                                                          &
1388           ONLY:  hyp, nr, pt, pt_init, prr, q, qc, qr, zu
1389
1390       USE cloud_parameters,                                                   &
1391           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
1392
1393       USE control_parameters,                                                 &
1394           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
1395                  intermediate_timestep_count, large_scale_forcing,            &
1396                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
1397                  pt_surface, rho_surface, surface_pressure
1398
1399       USE indices,                                                            &
1400           ONLY:  nzb, nzt
1401
1402       USE kinds
1403
1404       USE statistics,                                                         &
1405           ONLY:  weight_pres
1406
1407       IMPLICIT NONE
1408
1409       INTEGER(iwp) ::  i                 !<
1410       INTEGER(iwp) ::  j                 !<
1411       INTEGER(iwp) ::  k                 !<
1412
1413       REAL(wp)     ::  t_surface         !<
1414
1415       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
1416!
1417!--       Calculate:
1418!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1419!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1420!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1421          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1422          DO  k = nzb, nzt+1
1423             hyp(k)    = surface_pressure * 100.0_wp * &
1424                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1425             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1426             t_d_pt(k) = 1.0_wp / pt_d_t(k)
1427             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
1428          ENDDO
1429!
1430!--       Compute reference density
1431          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1432       ENDIF
1433
1434!
1435!--    Compute length of time step
1436       IF ( call_microphysics_at_all_substeps )  THEN
1437          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1438       ELSE
1439          dt_micro = dt_3d
1440       ENDIF
1441
1442!
1443!--    Use 1d arrays
1444       q_1d(:)  = q(:,j,i)
1445       pt_1d(:) = pt(:,j,i)
1446       qc_1d(:) = qc(:,j,i)
1447       nc_1d(:) = nc_const
1448       IF ( microphysics_seifert )  THEN
1449          qr_1d(:) = qr(:,j,i)
1450          nr_1d(:) = nr(:,j,i)
1451       ENDIF
1452
1453!
1454!--    Reset precipitation rate
1455       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1456
1457!
1458!--    Compute cloud physics
1459       IF( microphysics_kessler )  THEN
1460
1461          CALL autoconversion_kessler( i,j )
1462          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1463
1464       ELSEIF ( microphysics_seifert )  THEN
1465
1466          CALL adjust_cloud( i,j )
1467          CALL autoconversion( i,j )
1468          CALL accretion( i,j )
1469          CALL selfcollection_breakup( i,j )
1470          CALL evaporation_rain( i,j )
1471          CALL sedimentation_rain( i,j )
1472          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
1473
1474       ENDIF
1475
1476       CALL calc_precipitation_amount( i,j )
1477
1478!
1479!--    Store results on the 3d arrays
1480       q(:,j,i)  = q_1d(:)
1481       pt(:,j,i) = pt_1d(:)
1482       IF ( microphysics_seifert )  THEN
1483          qr(:,j,i) = qr_1d(:)
1484          nr(:,j,i) = nr_1d(:)
1485       ENDIF
1486
1487    END SUBROUTINE microphysics_control_ij
1488
1489!------------------------------------------------------------------------------!
1490! Description:
1491! ------------
1492!> Adjust number of raindrops to avoid nonlinear effects in
1493!> sedimentation and evaporation of rain drops due to too small or
1494!> too big weights of rain drops (Stevens and Seifert, 2008).
1495!> The same procedure is applied to cloud droplets if they are determined
1496!> prognostically. Call for grid point i,j
1497!------------------------------------------------------------------------------!
1498    SUBROUTINE adjust_cloud_ij( i, j )
1499
1500       USE cloud_parameters,                                                   &
1501           ONLY:  hyrho
1502
1503       USE indices,                                                            &
1504           ONLY:  nzb, nzt, wall_flags_0
1505
1506       USE kinds
1507
1508       IMPLICIT NONE
1509
1510       INTEGER(iwp) ::  i                 !<
1511       INTEGER(iwp) ::  j                 !<
1512       INTEGER(iwp) ::  k                 !<
1513
1514       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
1515
1516       DO  k = nzb+1, nzt
1517!
1518!--       Predetermine flag to mask topography
1519          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1520
1521          IF ( qr_1d(k) <= eps_sb )  THEN
1522             qr_1d(k) = 0.0_wp
1523             nr_1d(k) = 0.0_wp
1524          ELSE
1525!
1526!--          Adjust number of raindrops to avoid nonlinear effects in
1527!--          sedimentation and evaporation of rain drops due to too small or
1528!--          too big weights of rain drops (Stevens and Seifert, 2008).
1529             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
1530                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
1531             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
1532                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
1533             ENDIF
1534
1535          ENDIF
1536
1537       ENDDO
1538
1539    END SUBROUTINE adjust_cloud_ij
1540
1541
1542!------------------------------------------------------------------------------!
1543! Description:
1544! ------------
1545!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
1546!------------------------------------------------------------------------------!
1547    SUBROUTINE autoconversion_ij( i, j )
1548
1549       USE arrays_3d,                                                          &
1550           ONLY:  diss, dzu
1551
1552       USE cloud_parameters,                                                   &
1553           ONLY:  hyrho
1554
1555       USE control_parameters,                                                 &
1556           ONLY:  rho_surface
1557
1558       USE grid_variables,                                                     &
1559           ONLY:  dx, dy
1560
1561       USE indices,                                                            &
1562           ONLY:  nzb, nzt, wall_flags_0
1563
1564       USE kinds
1565
1566       IMPLICIT NONE
1567
1568       INTEGER(iwp) ::  i                 !<
1569       INTEGER(iwp) ::  j                 !<
1570       INTEGER(iwp) ::  k                 !<
1571
1572       REAL(wp)     ::  alpha_cc          !<
1573       REAL(wp)     ::  autocon           !<
1574       REAL(wp)     ::  dissipation       !<
1575       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
1576       REAL(wp)     ::  k_au              !<
1577       REAL(wp)     ::  l_mix             !<
1578       REAL(wp)     ::  nu_c              !<
1579       REAL(wp)     ::  phi_au            !<
1580       REAL(wp)     ::  r_cc              !<
1581       REAL(wp)     ::  rc                !<
1582       REAL(wp)     ::  re_lambda         !<
1583       REAL(wp)     ::  sigma_cc          !<
1584       REAL(wp)     ::  tau_cloud         !<
1585       REAL(wp)     ::  xc                !<
1586
1587       DO  k = nzb+1, nzt
1588!
1589!--       Predetermine flag to mask topography
1590          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1591
1592          IF ( qc_1d(k) > eps_sb )  THEN
1593
1594             k_au = k_cc / ( 20.0_wp * x0 )
1595!
1596!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1597!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1598             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
1599!
1600!--          Universal function for autoconversion process
1601!--          (Seifert and Beheng, 2006):
1602             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
1603!
1604!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
1605!--          (Use constant nu_c = 1.0_wp instead?)
1606             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
1607!
1608!--          Mean weight of cloud droplets:
1609             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
1610!
1611!--          Parameterized turbulence effects on autoconversion (Seifert,
1612!--          Nuijens and Stevens, 2010)
1613             IF ( collision_turbulence )  THEN
1614!
1615!--             Weight averaged radius of cloud droplets:
1616                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
1617
1618                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1619                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1620                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
1621!
1622!--             Mixing length (neglecting distance to ground and stratification)
1623                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
1624!
1625!--             Limit dissipation rate according to Seifert, Nuijens and
1626!--             Stevens (2010)
1627                dissipation = MIN( 0.06_wp, diss(k,j,i) )
1628!
1629!--             Compute Taylor-microscale Reynolds number:
1630                re_lambda = 6.0_wp / 11.0_wp *                                 &
1631                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
1632                            SQRT( 15.0_wp / kin_vis_air ) *                    &
1633                            dissipation**( 1.0_wp / 6.0_wp )
1634!
1635!--             The factor of 1.0E4 is needed to convert the dissipation rate
1636!--             from m2 s-3 to cm2 s-3.
1637                k_au = k_au * ( 1.0_wp +                                       &
1638                                dissipation * 1.0E4_wp *                       &
1639                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
1640                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
1641                                                  sigma_cc )**2                &
1642                                                ) + beta_cc                    &
1643                                )                                              &
1644                              )
1645             ENDIF
1646!
1647!--          Autoconversion rate (Seifert and Beheng, 2006):
1648             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
1649                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
1650                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
1651                       rho_surface
1652             autocon = MIN( autocon, qc_1d(k) / dt_micro )
1653
1654             qr_1d(k) = qr_1d(k) + autocon * dt_micro                 * flag
1655             qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
1656             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
1657
1658          ENDIF
1659
1660       ENDDO
1661
1662    END SUBROUTINE autoconversion_ij
1663
1664!------------------------------------------------------------------------------!
1665! Description:
1666! ------------
1667!> Autoconversion process (Kessler, 1969).
1668!------------------------------------------------------------------------------!
1669    SUBROUTINE autoconversion_kessler_ij( i, j )
1670
1671       USE arrays_3d,                                                          &
1672           ONLY:  dzw, prr
1673
1674       USE cloud_parameters,                                                   &
1675           ONLY:  l_d_cp, pt_d_t
1676
1677       USE indices,                                                            &
1678           ONLY:  nzb, nzb_max, nzt, wall_flags_0
1679
1680       USE kinds
1681
1682
1683       IMPLICIT NONE
1684
1685       INTEGER(iwp) ::  i      !<
1686       INTEGER(iwp) ::  j      !<
1687       INTEGER(iwp) ::  k      !<
1688       INTEGER(iwp) ::  k_wall !< topography top index
1689
1690       REAL(wp)    ::  dqdt_precip !<
1691       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
1692
1693!
1694!--    Determine vertical index of topography top
1695       k_wall = MAXLOC(                                                        &
1696                        MERGE( 1, 0,                                           &
1697                               BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )      &
1698                             ), DIM = 1                                        &
1699                      ) - 1
1700       DO  k = nzb+1, nzt
1701!
1702!--       Predetermine flag to mask topography
1703          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1704
1705          IF ( qc_1d(k) > ql_crit )  THEN
1706             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
1707          ELSE
1708             dqdt_precip = 0.0_wp
1709          ENDIF
1710
1711          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
1712          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro * flag
1713          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
1714
1715!
1716!--       Compute the rain rate (stored on surface grid point)
1717          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
1718
1719       ENDDO
1720
1721    END SUBROUTINE autoconversion_kessler_ij
1722
1723!------------------------------------------------------------------------------!
1724! Description:
1725! ------------
1726!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
1727!------------------------------------------------------------------------------!
1728    SUBROUTINE accretion_ij( i, j )
1729
1730       USE arrays_3d,                                                          &
1731           ONLY:  diss
1732
1733       USE cloud_parameters,                                                   &
1734           ONLY:  hyrho
1735
1736       USE control_parameters,                                                 &
1737           ONLY:  rho_surface
1738
1739       USE indices,                                                            &
1740           ONLY:  nzb, nzt, wall_flags_0
1741
1742       USE kinds
1743
1744       IMPLICIT NONE
1745
1746       INTEGER(iwp) ::  i                 !<
1747       INTEGER(iwp) ::  j                 !<
1748       INTEGER(iwp) ::  k                 !<
1749
1750       REAL(wp)     ::  accr              !<
1751       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
1752       REAL(wp)     ::  k_cr              !<
1753       REAL(wp)     ::  phi_ac            !<
1754       REAL(wp)     ::  tau_cloud         !<
1755
1756       DO  k = nzb+1, nzt
1757!
1758!--       Predetermine flag to mask topography
1759          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1760
1761          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
1762!
1763!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1764             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
1765!
1766!--          Universal function for accretion process
1767!--          (Seifert and Beheng, 2001):
1768             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1769!
1770!--          Parameterized turbulence effects on autoconversion (Seifert,
1771!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1772!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1773             IF ( collision_turbulence )  THEN
1774                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
1775                                 MIN( 600.0_wp,                                &
1776                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
1777                               )
1778             ELSE
1779                k_cr = k_cr0
1780             ENDIF
1781!
1782!--          Accretion rate (Seifert and Beheng, 2006):
1783             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
1784             accr = MIN( accr, qc_1d(k) / dt_micro )
1785
1786             qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
1787             qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
1788
1789          ENDIF
1790
1791       ENDDO
1792
1793    END SUBROUTINE accretion_ij
1794
1795
1796!------------------------------------------------------------------------------!
1797! Description:
1798! ------------
1799!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
1800!------------------------------------------------------------------------------!
1801    SUBROUTINE selfcollection_breakup_ij( i, j )
1802
1803       USE cloud_parameters,                                                   &
1804           ONLY:  hyrho
1805
1806       USE control_parameters,                                                 &
1807           ONLY:  rho_surface
1808
1809       USE indices,                                                            &
1810           ONLY:  nzb, nzt, wall_flags_0
1811
1812       USE kinds
1813
1814       IMPLICIT NONE
1815
1816       INTEGER(iwp) ::  i                 !<
1817       INTEGER(iwp) ::  j                 !<
1818       INTEGER(iwp) ::  k                 !<
1819
1820       REAL(wp)     ::  breakup           !<
1821       REAL(wp)     ::  dr                !<
1822       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
1823       REAL(wp)     ::  phi_br            !<
1824       REAL(wp)     ::  selfcoll          !<
1825
1826       DO  k = nzb+1, nzt
1827!
1828!--       Predetermine flag to mask topography
1829          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1830
1831          IF ( qr_1d(k) > eps_sb )  THEN
1832!
1833!--          Selfcollection rate (Seifert and Beheng, 2001):
1834             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
1835!
1836!--          Weight averaged diameter of rain drops:
1837             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1838!
1839!--          Collisional breakup rate (Seifert, 2008):
1840             IF ( dr >= 0.3E-3_wp )  THEN
1841                phi_br  = k_br * ( dr - 1.1E-3_wp )
1842                breakup = selfcoll * ( phi_br + 1.0_wp )
1843             ELSE
1844                breakup = 0.0_wp
1845             ENDIF
1846
1847             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1848             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
1849
1850          ENDIF
1851       ENDDO
1852
1853    END SUBROUTINE selfcollection_breakup_ij
1854
1855
1856!------------------------------------------------------------------------------!
1857! Description:
1858! ------------
1859!> Evaporation of precipitable water. Condensation is neglected for
1860!> precipitable water. Call for grid point i,j
1861!------------------------------------------------------------------------------!
1862    SUBROUTINE evaporation_rain_ij( i, j )
1863
1864       USE arrays_3d,                                                          &
1865           ONLY:  hyp
1866
1867       USE cloud_parameters,                                                   &
1868           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
1869
1870       USE constants,                                                          &
1871           ONLY:  pi
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)     ::  alpha             !<
1885       REAL(wp)     ::  dr                !<
1886       REAL(wp)     ::  e_s               !<
1887       REAL(wp)     ::  evap              !<
1888       REAL(wp)     ::  evap_nr           !<
1889       REAL(wp)     ::  f_vent            !<
1890       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
1891       REAL(wp)     ::  g_evap            !<
1892       REAL(wp)     ::  lambda_r          !<
1893       REAL(wp)     ::  mu_r              !<
1894       REAL(wp)     ::  mu_r_2            !<
1895       REAL(wp)     ::  mu_r_5d2          !<
1896       REAL(wp)     ::  nr_0              !<
1897       REAL(wp)     ::  q_s               !<
1898       REAL(wp)     ::  sat               !<
1899       REAL(wp)     ::  t_l               !<
1900       REAL(wp)     ::  temp              !<
1901       REAL(wp)     ::  xr                !<
1902
1903       DO  k = nzb+1, nzt
1904!
1905!--       Predetermine flag to mask topography
1906          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1907
1908          IF ( qr_1d(k) > eps_sb )  THEN
1909!
1910!--          Actual liquid water temperature:
1911             t_l = t_d_pt(k) * pt_1d(k)
1912!
1913!--          Saturation vapor pressure at t_l:
1914             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
1915                                    ( t_l - 35.86_wp )                         &
1916                                  )
1917!
1918!--          Computation of saturation humidity:
1919             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1920             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1921             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
1922!
1923!--          Supersaturation:
1924             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
1925!
1926!--          Evaporation needs only to be calculated in subsaturated regions
1927             IF ( sat < 0.0_wp )  THEN
1928!
1929!--             Actual temperature:
1930                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1931
1932                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
1933                                    ( thermal_conductivity_l * temp ) +        &
1934                                    r_v * temp / ( diff_coeff_l * e_s )        &
1935                                  )
1936!
1937!--             Mean weight of rain drops
1938                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
1939!
1940!--             Weight averaged diameter of rain drops:
1941                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1942!
1943!--             Compute ventilation factor and intercept parameter
1944!--             (Seifert and Beheng, 2006; Seifert, 2008):
1945                IF ( ventilation_effect )  THEN
1946!
1947!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1948!--                Stevens and Seifert, 2008):
1949                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1950!
1951!--                Slope parameter of gamma distribution (Seifert, 2008):
1952                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1953                                ( mu_r + 1.0_wp )                              &
1954                              )**( 1.0_wp / 3.0_wp ) / dr
1955
1956                   mu_r_2   = mu_r + 2.0_wp
1957                   mu_r_5d2 = mu_r + 2.5_wp
1958
1959                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
1960                            b_vent * schmidt_p_1d3 *                           &
1961                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
1962                            lambda_r**( -mu_r_5d2 ) *                          &
1963                            ( 1.0_wp -                                         &
1964                              0.5_wp * ( b_term / a_term ) *                   &
1965                              ( lambda_r / ( c_term + lambda_r )               &
1966                              )**mu_r_5d2 -                                    &
1967                              0.125_wp * ( b_term / a_term )**2 *              &
1968                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
1969                              )**mu_r_5d2 -                                    &
1970                              0.0625_wp * ( b_term / a_term )**3 *             &
1971                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
1972                              )**mu_r_5d2 -                                    &
1973                              0.0390625_wp * ( b_term / a_term )**4 *          &
1974                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
1975                              )**mu_r_5d2                                      &
1976                            )
1977
1978                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
1979                            gamm( mu_r + 1.0_wp )
1980                ELSE
1981                   f_vent = 1.0_wp
1982                   nr_0   = nr_1d(k) * dr
1983                ENDIF
1984!
1985!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
1986                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1987                evap    = MAX( evap, -qr_1d(k) / dt_micro )
1988                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
1989                               -nr_1d(k) / dt_micro )
1990
1991                qr_1d(k) = qr_1d(k) + evap    * dt_micro * flag
1992                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
1993
1994             ENDIF
1995          ENDIF
1996
1997       ENDDO
1998
1999    END SUBROUTINE evaporation_rain_ij
2000
2001
2002!------------------------------------------------------------------------------!
2003! Description:
2004! ------------
2005!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2006!> Call for grid point i,j
2007!------------------------------------------------------------------------------!
2008    SUBROUTINE sedimentation_cloud_ij( i, j )
2009
2010       USE arrays_3d,                                                          &
2011           ONLY:  ddzu, dzu, prr
2012
2013       USE cloud_parameters,                                                   &
2014           ONLY:  hyrho, l_d_cp, pt_d_t
2015
2016       USE control_parameters,                                                 &
2017           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
2018
2019       USE indices,                                                            &
2020           ONLY:  nzb, nzb, nzt, wall_flags_0
2021
2022       USE kinds
2023
2024       USE statistics,                                                         &
2025           ONLY:  weight_substep
2026
2027       IMPLICIT NONE
2028
2029       INTEGER(iwp) ::  i             !<
2030       INTEGER(iwp) ::  j             !<
2031       INTEGER(iwp) ::  k             !<
2032
2033       REAL(wp)                       ::  flag    !< flag to indicate first grid level above surface
2034       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
2035
2036       sed_qc(nzt+1) = 0.0_wp
2037
2038       DO  k = nzt, nzb+1, -1
2039!
2040!--       Predetermine flag to mask topography
2041          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2042
2043          IF ( qc_1d(k) > eps_sb )  THEN
2044             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
2045                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
2046          ELSE
2047             sed_qc(k) = 0.0_wp
2048          ENDIF
2049
2050          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
2051                                      dt_micro + sed_qc(k+1)                   &
2052                         ) * flag
2053
2054          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2055                                hyrho(k) * dt_micro * flag
2056          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2057                                hyrho(k) * dt_micro * flag
2058          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
2059                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
2060
2061!
2062!--       Compute the precipitation rate of cloud (fog) droplets
2063          IF ( call_microphysics_at_all_substeps )  THEN
2064             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
2065                              weight_substep(intermediate_timestep_count) * flag
2066          ELSE
2067             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2068          ENDIF
2069
2070       ENDDO
2071
2072    END SUBROUTINE sedimentation_cloud_ij
2073
2074
2075!------------------------------------------------------------------------------!
2076! Description:
2077! ------------
2078!> Computation of sedimentation flux. Implementation according to Stevens
2079!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
2080!------------------------------------------------------------------------------!
2081    SUBROUTINE sedimentation_rain_ij( i, j )
2082
2083       USE arrays_3d,                                                          &
2084           ONLY:  ddzu, dzu, prr
2085
2086       USE cloud_parameters,                                                   &
2087           ONLY:  hyrho, l_d_cp, pt_d_t
2088
2089       USE control_parameters,                                                 &
2090           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
2091
2092       USE indices,                                                            &
2093           ONLY:  nzb, nzb, nzt, wall_flags_0
2094
2095       USE kinds
2096
2097       USE statistics,                                                         &
2098           ONLY:  weight_substep
2099
2100       USE surface_mod,                                                        &
2101           ONLY :  bc_h
2102       
2103       IMPLICIT NONE
2104
2105       INTEGER(iwp) ::  i             !< running index x direction
2106       INTEGER(iwp) ::  j             !< running index y direction
2107       INTEGER(iwp) ::  k             !< running index z direction
2108       INTEGER(iwp) ::  k_run         !<
2109       INTEGER(iwp) ::  m             !< running index surface elements
2110       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2111       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2112
2113       REAL(wp)     ::  c_run                      !<
2114       REAL(wp)     ::  d_max                      !<
2115       REAL(wp)     ::  d_mean                     !<
2116       REAL(wp)     ::  d_min                      !<
2117       REAL(wp)     ::  dr                         !<
2118       REAL(wp)     ::  flux                       !<
2119       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
2120       REAL(wp)     ::  lambda_r                   !<
2121       REAL(wp)     ::  mu_r                       !<
2122       REAL(wp)     ::  z_run                      !<
2123
2124       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2125       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2126       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2127       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2128       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2129       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2130       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2131       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
2132
2133!
2134!--    Compute velocities
2135       DO  k = nzb+1, nzt
2136!
2137!--       Predetermine flag to mask topography
2138          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2139
2140          IF ( qr_1d(k) > eps_sb )  THEN
2141!
2142!--          Weight averaged diameter of rain drops:
2143             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
2144!
2145!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2146!--          Stevens and Seifert, 2008):
2147             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2148!
2149!--          Slope parameter of gamma distribution (Seifert, 2008):
2150             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
2151                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2152
2153             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2154                                         a_term - b_term * ( 1.0_wp +          &
2155                                            c_term / lambda_r )**( -1.0_wp *   &
2156                                            ( mu_r + 1.0_wp ) )                &
2157                                        )                                      &
2158                          ) * flag
2159             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2160                                         a_term - b_term * ( 1.0_wp +          &
2161                                            c_term / lambda_r )**( -1.0_wp *   &
2162                                            ( mu_r + 4.0_wp ) )                &
2163                                       )                                       &
2164                          ) * flag
2165          ELSE
2166             w_nr(k) = 0.0_wp
2167             w_qr(k) = 0.0_wp
2168          ENDIF
2169       ENDDO
2170!
2171!--    Adjust boundary values using surface data type.
2172!--    Upward facing non-natural
2173       surf_s = bc_h(0)%start_index(j,i)   
2174       surf_e = bc_h(0)%end_index(j,i)
2175       DO  m = surf_s, surf_e
2176          k         = bc_h(0)%k(m)
2177          w_nr(k-1) = w_nr(k)
2178          w_qr(k-1) = w_qr(k)
2179       ENDDO
2180!
2181!--    Downward facing non-natural
2182       surf_s = bc_h(1)%start_index(j,i)   
2183       surf_e = bc_h(1)%end_index(j,i)
2184       DO  m = surf_s, surf_e
2185          k         = bc_h(1)%k(m)
2186          w_nr(k+1) = w_nr(k)
2187          w_qr(k+1) = w_qr(k)
2188       ENDDO
2189!
2190!--    Neumann boundary condition at model top
2191       w_nr(nzt+1) = 0.0_wp
2192       w_qr(nzt+1) = 0.0_wp
2193!
2194!--    Compute Courant number
2195       DO  k = nzb+1, nzt
2196!
2197!--       Predetermine flag to mask topography
2198          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2199
2200          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
2201                    dt_micro * ddzu(k) * flag
2202          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
2203                    dt_micro * ddzu(k) * flag
2204       ENDDO
2205!
2206!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2207       IF ( limiter_sedimentation )  THEN
2208
2209          DO k = nzb+1, nzt
2210!
2211!--          Predetermine flag to mask topography
2212             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2213
2214             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
2215             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2216             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
2217
2218             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2219                                                        2.0_wp * d_max,        &
2220                                                        ABS( d_mean ) ) * flag
2221
2222             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
2223             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2224             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
2225
2226             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2227                                                        2.0_wp * d_max,        &
2228                                                        ABS( d_mean ) ) * flag
2229          ENDDO
2230
2231       ELSE
2232
2233          nr_slope = 0.0_wp
2234          qr_slope = 0.0_wp
2235
2236       ENDIF
2237
2238       sed_nr(nzt+1) = 0.0_wp
2239       sed_qr(nzt+1) = 0.0_wp
2240!
2241!--    Compute sedimentation flux
2242       DO  k = nzt, nzb+1, -1
2243!
2244!--       Predetermine flag to mask topography
2245          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2246!
2247!--       Sum up all rain drop number densities which contribute to the flux
2248!--       through k-1/2
2249          flux  = 0.0_wp
2250          z_run = 0.0_wp ! height above z(k)
2251          k_run = k
2252          c_run = MIN( 1.0_wp, c_nr(k) )
2253          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
2254             flux  = flux + hyrho(k_run) *                                     &
2255                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
2256                     0.5_wp ) * c_run * dzu(k_run) * flag
2257             z_run = z_run + dzu(k_run)            * flag
2258             k_run = k_run + 1                     * flag
2259             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
2260          ENDDO
2261!
2262!--       It is not allowed to sediment more rain drop number density than
2263!--       available
2264          flux = MIN( flux,                                                    &
2265                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
2266
2267          sed_nr(k) = flux / dt_micro * flag
2268          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
2269                                    hyrho(k) * dt_micro * flag
2270!
2271!--       Sum up all rain water content which contributes to the flux
2272!--       through k-1/2
2273          flux  = 0.0_wp
2274          z_run = 0.0_wp ! height above z(k)
2275          k_run = k
2276          c_run = MIN( 1.0_wp, c_qr(k) )
2277
2278          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
2279
2280             flux  = flux + hyrho(k_run) *                                     &
2281                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
2282                     0.5_wp ) * c_run * dzu(k_run) * flag
2283             z_run = z_run + dzu(k_run)            * flag
2284             k_run = k_run + 1                     * flag
2285             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
2286
2287          ENDDO
2288!
2289!--       It is not allowed to sediment more rain water content than available
2290          flux = MIN( flux,                                                    &
2291                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
2292
2293          sed_qr(k) = flux / dt_micro * flag
2294
2295          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2296                                hyrho(k) * dt_micro * flag
2297          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2298                                hyrho(k) * dt_micro * flag
2299          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
2300                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
2301!
2302!--       Compute the rain rate
2303          IF ( call_microphysics_at_all_substeps )  THEN
2304             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
2305                          * weight_substep(intermediate_timestep_count) * flag
2306          ELSE
2307             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
2308          ENDIF
2309
2310       ENDDO
2311
2312    END SUBROUTINE sedimentation_rain_ij
2313
2314
2315!------------------------------------------------------------------------------!
2316! Description:
2317! ------------
2318!> This subroutine computes the precipitation amount due to gravitational
2319!> settling of rain and cloud (fog) droplets
2320!------------------------------------------------------------------------------!
2321    SUBROUTINE calc_precipitation_amount_ij( i, j )
2322
2323       USE arrays_3d,                                                          &
2324           ONLY:  precipitation_amount, prr
2325
2326       USE cloud_parameters,                                                   &
2327           ONLY:  hyrho
2328
2329       USE control_parameters,                                                 &
2330           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
2331                  intermediate_timestep_count, intermediate_timestep_count_max,&
2332                  precipitation_amount_interval, time_do2d_xy
2333
2334       USE indices,                                                            &
2335           ONLY:  nzb, nzt, wall_flags_0
2336
2337       USE kinds
2338
2339       USE surface_mod,                                                        &
2340           ONLY :  bc_h
2341
2342       IMPLICIT NONE
2343
2344       INTEGER(iwp) ::  i             !< running index x direction
2345       INTEGER(iwp) ::  j             !< running index y direction
2346       INTEGER(iwp) ::  k             !< running index z direction
2347       INTEGER(iwp) ::  m             !< running index surface elements
2348       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2349       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2350
2351       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
2352            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
2353            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
2354       THEN
2355
2356          surf_s = bc_h(0)%start_index(j,i)   
2357          surf_e = bc_h(0)%end_index(j,i)
2358          DO  m = surf_s, surf_e
2359             k                         = bc_h(0)%k(m)
2360             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
2361                                               prr(k,j,i) * hyrho(k) * dt_3d
2362          ENDDO
2363
2364       ENDIF
2365
2366    END SUBROUTINE calc_precipitation_amount_ij
2367
2368!------------------------------------------------------------------------------!
2369! Description:
2370! ------------
2371!> This function computes the gamma function (Press et al., 1992).
2372!> The gamma function is needed for the calculation of the evaporation
2373!> of rain drops.
2374!------------------------------------------------------------------------------!
2375    FUNCTION gamm( xx )
2376
2377       USE kinds
2378
2379       IMPLICIT NONE
2380
2381       INTEGER(iwp) ::  j            !<
2382
2383       REAL(wp)     ::  gamm         !<
2384       REAL(wp)     ::  ser          !<
2385       REAL(wp)     ::  tmp          !<
2386       REAL(wp)     ::  x_gamm       !<
2387       REAL(wp)     ::  xx           !<
2388       REAL(wp)     ::  y_gamm       !<
2389
2390
2391       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
2392       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
2393                                           -86.50532032941677_wp,      &
2394                                            24.01409824083091_wp,      &
2395                                            -1.231739572450155_wp,     &
2396                                             0.1208650973866179E-2_wp, &
2397                                            -0.5395239384953E-5_wp /)     !<
2398
2399       x_gamm = xx
2400       y_gamm = x_gamm
2401       tmp = x_gamm + 5.5_wp
2402       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
2403       ser = 1.000000000190015_wp
2404
2405       DO  j = 1, 6
2406          y_gamm = y_gamm + 1.0_wp
2407          ser    = ser + cof( j ) / y_gamm
2408       ENDDO
2409
2410!
2411!--    Until this point the algorithm computes the logarithm of the gamma
2412!--    function. Hence, the exponential function is used.
2413!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
2414       gamm = EXP( tmp ) * stp * ser / x_gamm
2415
2416       RETURN
2417
2418    END FUNCTION gamm
2419
2420 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.