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

Last change on this file since 2345 was 2318, checked in by suehring, 7 years ago

get topograpyhy top index via function call

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