source: palm/tags/release-5.0/SOURCE/microphysics_mod.f90 @ 2704

Last change on this file since 2704 was 2701, checked in by suehring, 6 years ago

changes from last commit documented

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