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

Last change on this file since 2307 was 2292, checked in by schwenkel, 7 years ago

implementation of new bulk microphysics scheme

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