source: palm/trunk/SOURCE/bulk_cloud_model_mod.f90 @ 3860

Last change on this file since 3860 was 3786, checked in by raasch, 5 years ago

unused variables removed, interoperable C datatypes introduced in particle type to avoid compiler warnings

  • Property svn:keywords set to Id
File size: 160.6 KB
Line 
1!> @file bulk_cloud_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: bulk_cloud_model_mod.f90 3786 2019-03-06 16:58:03Z raasch $
27! unsed variables removed
28!
29! 3767 2019-02-27 08:18:02Z raasch
30! unused variable for file index removed from rrd-subroutines parameter list
31!
32! 3724 2019-02-06 16:28:23Z kanani
33! Correct double-used log_point_s unit
34!
35! 3700 2019-01-26 17:03:42Z knoop
36! nopointer option removed
37!
38! 3622 2018-12-12 09:52:53Z schwenkel
39! Important bugfix in case of restart runs.
40!
41! 3507 2018-11-08 14:28:46Z schwenkel
42! Minor bugfixes for bcm cache version and initializing prr
43!
44! 3452 2018-10-30 13:13:34Z schwenkel
45! Bugfix for profiles output
46!
47! 3445 2018-10-29 12:23:02Z schwenkel
48! Minor bugfix and use of subroutine for supersaturation calculation in case
49! of cache version
50!
51! 3383 2018-10-19 14:22:58Z knoop
52! Modularization of all bulk cloud physics code components
53!
54! unused variables removed
55!
56! 3026 2018-05-22 10:30:53Z schwenkel
57! Changed the name specific humidity to mixing ratio, since we are computing
58! mixing ratios.
59!
60! 2718 2018-01-02 08:49:38Z maronga
61! Corrected "Former revisions" section
62!
63! 2701 2017-12-15 15:40:50Z suehring
64! Changes from last commit documented
65!
66! 2698 2017-12-14 18:46:24Z suehring
67! Bugfix in get_topography_top_index
68!
69! 2696 2017-12-14 17:12:51Z kanani
70! Change in file header (GPL part)
71!
72! 2608 2017-11-13 14:04:26Z schwenkel
73! Calculation of supersaturation in external module (diagnostic_quantities_mod).
74! Change: correct calculation of saturation specific humidity to saturation
75! mixing ratio (the factor of 0.378 vanishes).
76!
77! 2522 2017-10-05 14:20:37Z schwenkel
78! Minor bugfix
79!
80! 2375 2017-08-29 14:10:28Z schwenkel
81! Improved aerosol initilization and some minor bugfixes
82! for droplet sedimenation
83!
84! 2318 2017-07-20 17:27:44Z suehring
85! Get topography top index via Function call
86!
87! 2317 2017-07-20 17:27:19Z suehring
88! s1 changed to log_sigma
89!
90! 2292 2017-06-20 09:51:42Z schwenkel
91! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
92! includes two more prognostic equations for cloud drop concentration (nc)
93! and cloud water content (qc).
94! - The process of activation is parameterized with a simple Twomey
95!   activion scheme or with considering solution and curvature
96!   effects (Khvorostyanov and Curry ,2006).
97! - The saturation adjustment scheme is replaced by the parameterization
98!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
99! - All other microphysical processes of Seifert and Beheng are used.
100!   Additionally, in those processes the reduction of cloud number concentration
101!   is considered.
102!
103! 2233 2017-05-30 18:08:54Z suehring
104!
105! 2232 2017-05-30 17:47:52Z suehring
106! Adjustments to new topography and surface concept
107!
108! 2155 2017-02-21 09:57:40Z hoffmann
109! Bugfix in the calculation of microphysical quantities on ghost points.
110!
111! 2031 2016-10-21 15:11:58Z knoop
112! renamed variable rho to rho_ocean
113!
114! 2000 2016-08-20 18:09:15Z knoop
115! Forced header and separation lines into 80 columns
116!
117! 1850 2016-04-08 13:29:27Z maronga
118! Module renamed
119! Adapted for modularization of microphysics.
120!
121! 1845 2016-04-08 08:29:13Z raasch
122! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
123! point (instead of one point above surface)
124!
125! 1831 2016-04-07 13:15:51Z hoffmann
126! turbulence renamed collision_turbulence,
127! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
128! avaialble for microphysics_kessler.
129!
130! 1822 2016-04-07 07:49:42Z hoffmann
131! Unused variables removed.
132! Kessler scheme integrated.
133!
134! 1691 2015-10-26 16:17:44Z maronga
135! Added new routine calc_precipitation_amount. The routine now allows to account
136! for precipitation due to sedimenation of cloud (fog) droplets
137!
138! 1682 2015-10-07 23:56:08Z knoop
139! Code annotations made doxygen readable
140!
141! 1646 2015-09-02 16:00:10Z hoffmann
142! Bugfix: Wrong computation of d_mean.
143!
144! 1361 2014-04-16 15:17:48Z hoffmann
145! Bugfix in sedimentation_rain: Index corrected.
146! Vectorized version of adjust_cloud added.
147! Little reformatting of the code.
148!
149! 1353 2014-04-08 15:21:23Z heinze
150! REAL constants provided with KIND-attribute
151!
152! 1346 2014-03-27 13:18:20Z heinze
153! Bugfix: REAL constants provided with KIND-attribute especially in call of
154! intrinsic function like MAX, MIN, SIGN
155!
156! 1334 2014-03-25 12:21:40Z heinze
157! Bugfix: REAL constants provided with KIND-attribute
158!
159! 1322 2014-03-20 16:38:49Z raasch
160! REAL constants defined as wp-kind
161!
162! 1320 2014-03-20 08:40:49Z raasch
163! ONLY-attribute added to USE-statements,
164! kind-parameters added to all INTEGER and REAL declaration statements,
165! kinds are defined in new module kinds,
166! comment fields (!:) to be used for variable explanations added to
167! all variable declaration statements
168!
169! 1241 2013-10-30 11:36:58Z heinze
170! hyp and rho_ocean have to be calculated at each time step if data from external
171! file LSF_DATA are used
172!
173! 1115 2013-03-26 18:16:16Z hoffmann
174! microphyical tendencies are calculated in bcm_actions in an optimized
175! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
176!
177! 1106 2013-03-04 05:31:38Z raasch
178! small changes in code formatting
179!
180! 1092 2013-02-02 11:24:22Z raasch
181! unused variables removed
182! file put under GPL
183!
184! 1065 2012-11-22 17:42:36Z hoffmann
185! Sedimentation process implemented according to Stevens and Seifert (2008).
186! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
187! and Stevens, 2010).
188!
189! 1053 2012-11-13 17:11:03Z hoffmann
190! initial revision
191!
192! Description:
193! ------------
194!> Calculate bulk cloud microphysics.
195!------------------------------------------------------------------------------!
196 MODULE bulk_cloud_model_mod
197
198    USE arrays_3d,                                                             &
199        ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
200               nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p,         &
201               precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
202               qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
203               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
204
205    USE averaging,                                                             &
206        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
207
208    USE basic_constants_and_equations_mod,                                     &
209        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
210               molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
211               exner_function, exner_function_invers, ideal_gas_law_rho,       &
212               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
213
214    USE control_parameters,                                                    &
215        ONLY:  dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
216               intermediate_timestep_count_max, large_scale_forcing,           &
217               lsf_surf, pt_surface, rho_surface, surface_pressure,            &
218               time_do2d_xy, message_string, initializing_actions
219
220    USE cpulog,                                                                &
221        ONLY:  cpu_log, log_point_s
222
223    USE grid_variables,                                                        &
224        ONLY:  dx, dy
225
226    USE indices,                                                               &
227        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,           &
228               wall_flags_0
229
230    USE  kinds
231
232    USE statistics,                                                            &
233        ONLY:  weight_pres, weight_substep
234
235    USE surface_mod,                                                           &
236        ONLY :  bc_h, get_topography_top_index_ji, surf_bulk_cloud_model,      &
237                surf_microphysics_morrison, surf_microphysics_seifert
238
239    IMPLICIT NONE
240
241    CHARACTER (LEN=20)   ::  aerosol_bulk = 'nacl'                        !< namelist parameter
242    CHARACTER (LEN=20)   ::  cloud_scheme = 'saturation_adjust'           !< namelist parameter
243
244    LOGICAL ::  aerosol_nacl =.TRUE.                             !< nacl aerosol for bulk scheme
245    LOGICAL ::  aerosol_c3h4o4 =.FALSE.                          !< malonic acid aerosol for bulk scheme
246    LOGICAL ::  aerosol_nh4no3 =.FALSE.                          !< malonic acid aerosol for bulk scheme
247
248    LOGICAL ::  bulk_cloud_model = .FALSE.                       !< namelist parameter
249
250    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
251    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
252    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
253    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
254    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
255
256    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
257    LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
258    LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
259    LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
260    LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
261    LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
262
263    REAL(wp) ::  precipitation_amount_interval = 9999999.9_wp    !< namelist parameter
264
265    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
266    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
267    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
268    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
269    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
270    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
271    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
272    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
273    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
274    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
275    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
276    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
277    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
278    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
279    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
280    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
281    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
282    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
283    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
284    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
285    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
286    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
287    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
288    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
289    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
290    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
291    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
292    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
293    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
294    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
295    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
296    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
297    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
298    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
299    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
300    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
301    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
302
303    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
304    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
305    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
306    REAL(wp) ::  dt_micro                        !< microphysics time step
307    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
308    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
309    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
310    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
311    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
312    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
313
314    REAL(wp) ::  e_s     !< saturation water vapor pressure
315    REAL(wp) ::  q_s     !< saturation mixing ratio
316    REAL(wp) ::  sat     !< supersaturation
317    REAL(wp) ::  t_l     !< actual temperature
318
319    SAVE
320
321    PRIVATE
322
323    PUBLIC bcm_parin, &
324           bcm_check_parameters, &
325           bcm_check_data_output, &
326           bcm_check_data_output_pr, &
327           bcm_header, &
328           bcm_init_arrays, &
329           bcm_init, &
330           bcm_3d_data_averaging, &
331           bcm_data_output_2d, &
332           bcm_data_output_3d, &
333           bcm_swap_timelevel, &
334           bcm_rrd_global, &
335           bcm_rrd_local, &
336           bcm_wrd_global, &
337           bcm_wrd_local, &
338           bcm_actions, &
339           calc_liquid_water_content
340
341    PUBLIC call_microphysics_at_all_substeps, &
342           cloud_water_sedimentation, &
343           bulk_cloud_model, &
344           cloud_scheme, &
345           collision_turbulence, &
346           dt_precipitation, &
347           microphysics_morrison, &
348           microphysics_sat_adjust, &
349           microphysics_seifert, &
350           na_init, &
351           nc_const, &
352           precipitation, &
353           sigma_gc
354
355
356    INTERFACE bcm_parin
357       MODULE PROCEDURE bcm_parin
358    END INTERFACE bcm_parin
359
360    INTERFACE bcm_check_parameters
361       MODULE PROCEDURE bcm_check_parameters
362    END INTERFACE bcm_check_parameters
363
364    INTERFACE bcm_check_data_output
365       MODULE PROCEDURE bcm_check_data_output
366    END INTERFACE bcm_check_data_output
367
368    INTERFACE bcm_check_data_output_pr
369       MODULE PROCEDURE bcm_check_data_output_pr
370    END INTERFACE bcm_check_data_output_pr
371
372    INTERFACE bcm_header
373       MODULE PROCEDURE bcm_header
374    END INTERFACE bcm_header
375
376    INTERFACE bcm_init_arrays
377       MODULE PROCEDURE bcm_init_arrays
378    END INTERFACE bcm_init_arrays
379
380    INTERFACE bcm_init
381       MODULE PROCEDURE bcm_init
382    END INTERFACE bcm_init
383
384    INTERFACE bcm_3d_data_averaging
385       MODULE PROCEDURE bcm_3d_data_averaging
386    END INTERFACE bcm_3d_data_averaging
387
388    INTERFACE bcm_data_output_2d
389       MODULE PROCEDURE bcm_data_output_2d
390    END INTERFACE bcm_data_output_2d
391
392    INTERFACE bcm_data_output_3d
393       MODULE PROCEDURE bcm_data_output_3d
394    END INTERFACE bcm_data_output_3d
395
396    INTERFACE bcm_swap_timelevel
397       MODULE PROCEDURE bcm_swap_timelevel
398    END INTERFACE bcm_swap_timelevel
399
400    INTERFACE bcm_rrd_global
401       MODULE PROCEDURE bcm_rrd_global
402    END INTERFACE bcm_rrd_global
403
404    INTERFACE bcm_rrd_local
405       MODULE PROCEDURE bcm_rrd_local
406    END INTERFACE bcm_rrd_local
407
408    INTERFACE bcm_wrd_global
409       MODULE PROCEDURE bcm_wrd_global
410    END INTERFACE bcm_wrd_global
411
412    INTERFACE bcm_wrd_local
413       MODULE PROCEDURE bcm_wrd_local
414    END INTERFACE bcm_wrd_local
415
416    INTERFACE bcm_actions
417       MODULE PROCEDURE bcm_actions
418       MODULE PROCEDURE bcm_actions_ij
419    END INTERFACE bcm_actions
420
421    INTERFACE adjust_cloud
422       MODULE PROCEDURE adjust_cloud
423       MODULE PROCEDURE adjust_cloud_ij
424    END INTERFACE adjust_cloud
425
426    INTERFACE activation
427       MODULE PROCEDURE activation
428       MODULE PROCEDURE activation_ij
429    END INTERFACE activation
430
431    INTERFACE condensation
432       MODULE PROCEDURE condensation
433       MODULE PROCEDURE condensation_ij
434    END INTERFACE condensation
435
436    INTERFACE autoconversion
437       MODULE PROCEDURE autoconversion
438       MODULE PROCEDURE autoconversion_ij
439    END INTERFACE autoconversion
440
441    INTERFACE autoconversion_kessler
442       MODULE PROCEDURE autoconversion_kessler
443       MODULE PROCEDURE autoconversion_kessler_ij
444    END INTERFACE autoconversion_kessler
445
446    INTERFACE accretion
447       MODULE PROCEDURE accretion
448       MODULE PROCEDURE accretion_ij
449    END INTERFACE accretion
450
451    INTERFACE selfcollection_breakup
452       MODULE PROCEDURE selfcollection_breakup
453       MODULE PROCEDURE selfcollection_breakup_ij
454    END INTERFACE selfcollection_breakup
455
456    INTERFACE evaporation_rain
457       MODULE PROCEDURE evaporation_rain
458       MODULE PROCEDURE evaporation_rain_ij
459    END INTERFACE evaporation_rain
460
461    INTERFACE sedimentation_cloud
462       MODULE PROCEDURE sedimentation_cloud
463       MODULE PROCEDURE sedimentation_cloud_ij
464    END INTERFACE sedimentation_cloud
465
466    INTERFACE sedimentation_rain
467       MODULE PROCEDURE sedimentation_rain
468       MODULE PROCEDURE sedimentation_rain_ij
469    END INTERFACE sedimentation_rain
470
471    INTERFACE calc_precipitation_amount
472       MODULE PROCEDURE calc_precipitation_amount
473       MODULE PROCEDURE calc_precipitation_amount_ij
474    END INTERFACE calc_precipitation_amount
475
476    INTERFACE supersaturation
477       MODULE PROCEDURE supersaturation
478    END INTERFACE supersaturation
479
480    INTERFACE calc_liquid_water_content
481       MODULE PROCEDURE calc_liquid_water_content
482    END INTERFACE calc_liquid_water_content
483
484 CONTAINS
485
486
487!------------------------------------------------------------------------------!
488! Description:
489! ------------
490!> Parin for &bulk_cloud_parameters for the bulk cloud module
491!------------------------------------------------------------------------------!
492    SUBROUTINE bcm_parin
493
494
495       IMPLICIT NONE
496
497       CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
498
499       NAMELIST /bulk_cloud_parameters/  &
500          aerosol_bulk, &
501          c_sedimentation, &
502          call_microphysics_at_all_substeps, &
503          bulk_cloud_model, &
504          cloud_scheme, &
505          cloud_water_sedimentation, &
506          collision_turbulence, &
507          curvature_solution_effects_bulk, &
508          dry_aerosol_radius, &
509          limiter_sedimentation, &
510          na_init, &
511          nc_const, &
512          precipitation, &
513          precipitation_amount_interval, &
514          sigma_bulk, &
515          ventilation_effect
516
517       line = ' '
518!
519!--    Try to find bulk cloud module namelist
520       REWIND ( 11 )
521       line = ' '
522       DO   WHILE ( INDEX( line, '&bulk_cloud_parameters' ) == 0 )
523          READ ( 11, '(A)', END=10 )  line
524       ENDDO
525       BACKSPACE ( 11 )
526!
527!--    Read user-defined namelist
528       READ ( 11, bulk_cloud_parameters )
529!
530!--    Set flag that indicates that the bulk cloud module is switched on
531       !bulk_cloud_model = .TRUE.
532
53310     CONTINUE
534
535
536    END SUBROUTINE bcm_parin
537
538
539!------------------------------------------------------------------------------!
540! Description:
541! ------------
542!> Check parameters routine for bulk cloud module
543!------------------------------------------------------------------------------!
544    SUBROUTINE bcm_check_parameters
545
546
547       IMPLICIT NONE
548!
549!--    Check cloud scheme
550       IF ( cloud_scheme == 'saturation_adjust' )  THEN
551          microphysics_sat_adjust = .TRUE.
552          microphysics_seifert    = .FALSE.
553          microphysics_kessler    = .FALSE.
554          precipitation           = .FALSE.
555       ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
556          microphysics_sat_adjust = .FALSE.
557          microphysics_seifert    = .TRUE.
558          microphysics_kessler    = .FALSE.
559          microphysics_morrison  = .FALSE.
560          precipitation           = .TRUE.
561       ELSEIF ( cloud_scheme == 'kessler' )  THEN
562          microphysics_sat_adjust = .FALSE.
563          microphysics_seifert    = .FALSE.
564          microphysics_kessler    = .TRUE.
565          microphysics_morrison   = .FALSE.
566          precipitation           = .TRUE.
567       ELSEIF ( cloud_scheme == 'morrison' )  THEN
568          microphysics_sat_adjust = .FALSE.
569          microphysics_seifert    = .TRUE.
570          microphysics_kessler    = .FALSE.
571          microphysics_morrison   = .TRUE.
572          precipitation           = .TRUE.
573       ELSE
574          message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
575                           TRIM( cloud_scheme ) // '"'
576          CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
577       ENDIF
578
579
580
581!
582!--    Set the default value for the integration interval of precipitation amount
583       IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
584          IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
585             precipitation_amount_interval = dt_do2d_xy
586          ELSE
587             IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
588                WRITE( message_string, * )  'precipitation_amount_interval = ',   &
589                    precipitation_amount_interval, ' must not be larger than ',   &
590                    'dt_do2d_xy = ', dt_do2d_xy
591                CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
592             ENDIF
593          ENDIF
594       ENDIF
595
596       ! TODO: find better sollution for circular dependency problem
597       surf_bulk_cloud_model = bulk_cloud_model
598       surf_microphysics_morrison = microphysics_morrison
599       surf_microphysics_seifert = microphysics_seifert
600
601!
602!--    Check aerosol
603       IF ( aerosol_bulk == 'nacl' )  THEN
604          aerosol_nacl   = .TRUE.
605          aerosol_c3h4o4 = .FALSE.
606          aerosol_nh4no3 = .FALSE.
607       ELSEIF ( aerosol_bulk == 'c3h4o4' )  THEN
608          aerosol_nacl   = .FALSE.
609          aerosol_c3h4o4 = .TRUE.
610          aerosol_nh4no3 = .FALSE.
611       ELSEIF ( aerosol_bulk == 'nh4no3' )  THEN
612          aerosol_nacl   = .FALSE.
613          aerosol_c3h4o4 = .FALSE.
614          aerosol_nh4no3 = .TRUE.
615       ELSE
616          message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"'
617          CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 )
618       ENDIF
619
620
621    END SUBROUTINE bcm_check_parameters
622
623!------------------------------------------------------------------------------!
624! Description:
625! ------------
626!> Check data output for bulk cloud module
627!------------------------------------------------------------------------------!
628    SUBROUTINE bcm_check_data_output( var, unit )
629
630       IMPLICIT NONE
631
632       CHARACTER (LEN=*) ::  unit  !<
633       CHARACTER (LEN=*) ::  var   !<
634
635       SELECT CASE ( TRIM( var ) )
636
637          CASE ( 'nc' )
638             IF ( .NOT.  microphysics_morrison )  THEN
639                message_string = 'output of "' // TRIM( var ) // '" ' //       &
640                                 'requires ' //                                &
641                                 'cloud_scheme = "morrison"'
642                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
643             ENDIF
644             unit = '1/m3'
645
646          CASE ( 'nr' )
647             IF ( .NOT.  microphysics_seifert )  THEN
648                message_string = 'output of "' // TRIM( var ) // '" ' //       &
649                                 'requires ' //                                &
650                                 'cloud_scheme = "seifert_beheng"'
651                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
652             ENDIF
653             unit = '1/m3'
654
655          CASE ( 'prr' )
656             IF ( microphysics_sat_adjust )  THEN
657                message_string = 'output of "' // TRIM( var ) // '" ' //       &
658                                 'is not available for ' //                    &
659                                 'cloud_scheme = "saturation_adjust"'
660                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
661             ENDIF
662             unit = 'kg/kg m/s'
663
664          CASE ( 'qc' )
665             unit = 'kg/kg'
666
667          CASE ( 'qr' )
668             IF ( .NOT.  microphysics_seifert ) THEN
669                message_string = 'output of "' // TRIM( var ) // '" ' //       &
670                                 'requires ' //                                &
671                                 'cloud_scheme = "seifert_beheng"'
672                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
673             ENDIF
674             unit = 'kg/kg'
675
676          CASE ( 'pra*' )
677             IF ( .NOT. microphysics_kessler  .AND.                            &
678                  .NOT. microphysics_seifert )  THEN
679                message_string = 'output of "' // TRIM( var ) // '" ' //       &
680                                 'requires ' //                                &
681                                 'cloud_scheme = "kessler" or "seifert_beheng"'
682                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
683             ENDIF
684! TODO: find sollution (maybe connected to flow_statistics redesign?)
685!              IF ( j == 1 )  THEN
686!                 message_string = 'temporal averaging of precipitation ' //     &
687!                           'amount "' // TRIM( var ) // '" is not possible'
688!                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
689!              ENDIF
690             unit = 'mm'
691
692          CASE ( 'prr*' )
693             IF ( .NOT. microphysics_kessler  .AND.                            &
694                  .NOT. microphysics_seifert )  THEN
695                message_string = 'output of "' // TRIM( var ) // '"' //        &
696                         ' requires' //                                        &
697                         ' cloud_scheme = "kessler" or "seifert_beheng"'
698                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
699             ENDIF
700             unit = 'mm/s'
701
702          CASE DEFAULT
703             unit = 'illegal'
704
705       END SELECT
706
707
708    END SUBROUTINE bcm_check_data_output
709
710
711!------------------------------------------------------------------------------!
712! Description:
713! ------------
714!> Check data output of profiles for bulk cloud module
715!------------------------------------------------------------------------------!
716    SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit )
717
718       USE arrays_3d,                                                          &
719           ONLY: zu
720
721       USE control_parameters,                                                 &
722           ONLY: data_output_pr
723
724       USE profil_parameter,                                                   &
725           ONLY: dopr_index
726
727       USE statistics,                                                         &
728           ONLY: hom, statistic_regions
729
730       IMPLICIT NONE
731
732       CHARACTER (LEN=*) ::  unit      !<
733       CHARACTER (LEN=*) ::  variable  !<
734       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
735
736       INTEGER(iwp) ::  var_count      !<
737       INTEGER(iwp) ::  pr_index       !<
738
739       SELECT CASE ( TRIM( variable ) )
740
741! TODO: make index generic: pr_index = pr_palm+1
742
743          CASE ( 'nc' )
744             IF ( .NOT.  microphysics_morrison )  THEN
745                message_string = 'data_output_pr = ' //                        &
746                                 TRIM( data_output_pr(var_count) ) //          &
747                                 ' is not implemented for' //                  &
748                                 ' cloud_scheme /= morrison'
749                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
750             ENDIF
751             pr_index = 123
752             dopr_index(var_count) = pr_index
753             dopr_unit     = '1/m3'
754             unit = dopr_unit
755             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
756
757          CASE ( 'nr' )
758             IF ( .NOT.  microphysics_seifert )  THEN
759                message_string = 'data_output_pr = ' //                        &
760                                 TRIM( data_output_pr(var_count) ) //          &
761                                 ' is not implemented for' //                  &
762                                 ' cloud_scheme /= seifert_beheng'
763                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
764             ENDIF
765             pr_index = 73
766             dopr_index(var_count) = pr_index
767             dopr_unit     = '1/m3'
768             unit = dopr_unit
769             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
770
771          CASE ( 'prr' )
772             IF ( microphysics_sat_adjust )  THEN
773                message_string = 'data_output_pr = ' //                        &
774                                 TRIM( data_output_pr(var_count) ) //          &
775                                 ' is not available for' //                    &
776                                 ' cloud_scheme = saturation_adjust'
777                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
778             ENDIF
779             pr_index = 76
780             dopr_index(var_count) = pr_index
781             dopr_unit     = 'kg/kg m/s'
782             unit = dopr_unit
783             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
784          CASE ( 'qc' )
785             pr_index = 75
786             dopr_index(var_count) = pr_index
787             dopr_unit     = 'kg/kg'
788             unit = dopr_unit
789             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
790
791          CASE ( 'qr' )
792             IF ( .NOT.  microphysics_seifert )  THEN
793                message_string = 'data_output_pr = ' //                        &
794                                 TRIM( data_output_pr(var_count) ) //          &
795                                 ' is not implemented for' //                  &
796                                 ' cloud_scheme /= seifert_beheng'
797                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
798             ENDIF
799             pr_index = 74
800             dopr_index(var_count) = pr_index
801             dopr_unit     = 'kg/kg'
802             unit = dopr_unit
803             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
804
805          CASE DEFAULT
806             unit = 'illegal'
807
808       END SELECT
809
810    END SUBROUTINE bcm_check_data_output_pr
811
812
813!------------------------------------------------------------------------------!
814! Description:
815! ------------
816!> Allocate bulk cloud module arrays and define pointers
817!------------------------------------------------------------------------------!
818    SUBROUTINE bcm_init_arrays
819
820       USE indices,                                                            &
821           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
822
823
824       IMPLICIT NONE
825
826!
827!--    Liquid water content
828       ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
829
830!
831!--    3D-cloud water content
832       IF ( .NOT. microphysics_morrison )  THEN
833          ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
834       ENDIF
835!
836!--    Precipitation amount and rate (only needed if output is switched)
837       ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) )
838
839!
840!--    3d-precipitation rate
841       ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
842
843       IF ( microphysics_morrison )  THEN
844!
845!--       3D-cloud drop water content, cloud drop concentration arrays
846          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
847                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
848                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
849                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
850                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
851                    qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
852       ENDIF
853
854       IF ( microphysics_seifert )  THEN
855!
856!--       3D-rain water content, rain drop concentration arrays
857          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
858                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
859                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
860                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
861                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
862                    qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
863       ENDIF
864
865!
866!--    Initial assignment of the pointers
867       ql => ql_1
868       IF ( .NOT. microphysics_morrison )  THEN
869          qc => qc_1
870       ENDIF
871       IF ( microphysics_morrison )  THEN
872          qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
873          nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
874       ENDIF
875       IF ( microphysics_seifert )  THEN
876          qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
877          nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
878       ENDIF
879
880
881    END SUBROUTINE bcm_init_arrays
882
883
884!------------------------------------------------------------------------------!
885! Description:
886! ------------
887!> Initialization of the bulk cloud module
888!------------------------------------------------------------------------------!
889    SUBROUTINE bcm_init
890
891       IMPLICIT NONE
892
893       INTEGER(iwp) ::  i !<
894       INTEGER(iwp) ::  j !<
895
896       CALL location_message( 'initializing bulk cloud module', .FALSE. )
897
898       IF ( bulk_cloud_model )  THEN
899          IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
900!
901!--          Initialize the remaining quantities
902             IF ( microphysics_morrison )  THEN
903                DO  i = nxlg, nxrg
904                   DO  j = nysg, nyng
905                      qc(:,j,i) = 0.0_wp
906                      nc(:,j,i) = 0.0_wp
907                   ENDDO
908                ENDDO
909             ENDIF
910
911             IF ( microphysics_seifert )  THEN
912                DO  i = nxlg, nxrg
913                   DO  j = nysg, nyng
914                      qr(:,j,i) = 0.0_wp
915                      nr(:,j,i) = 0.0_wp
916                   ENDDO
917                ENDDO
918             ENDIF
919!
920!--          Liquid water content and precipitation amount
921!--          are zero at beginning of the simulation
922             ql = 0.0_wp
923             qc = 0.0_wp
924             precipitation_amount = 0.0_wp
925             prr = 0.0_wp
926!
927!--          Initialize old and new time levels.
928             IF ( microphysics_morrison )  THEN
929                tqc_m = 0.0_wp
930                tnc_m = 0.0_wp
931                qc_p  = qc
932                nc_p  = nc
933             ENDIF
934             IF ( microphysics_seifert )  THEN
935                tqr_m = 0.0_wp
936                tnr_m = 0.0_wp
937                qr_p  = qr
938                nr_p  = nr
939             ENDIF
940          ENDIF ! Only if not read_restart_data
941!
942!--       constant for the sedimentation of cloud water (2-moment cloud physics)
943          sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
944                             )**( 2.0_wp / 3.0_wp ) *                             &
945                         EXP( 5.0_wp * LOG( sigma_gc )**2 )
946
947!
948!--       Calculate timestep according to precipitation
949          IF ( microphysics_seifert )  THEN
950             dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
951                                w_precipitation
952          ENDIF
953
954!
955!--       Set constants for certain aerosol type
956          IF ( microphysics_morrison )  THEN
957             IF ( aerosol_nacl ) THEN
958                molecular_weight_of_solute = 0.05844_wp
959                rho_s                      = 2165.0_wp
960                vanthoff                   = 2.0_wp
961             ELSEIF ( aerosol_c3h4o4 ) THEN
962                molecular_weight_of_solute = 0.10406_wp
963                rho_s                      = 1600.0_wp
964                vanthoff                   = 1.37_wp
965             ELSEIF ( aerosol_nh4no3 ) THEN
966                molecular_weight_of_solute = 0.08004_wp
967                rho_s                      = 1720.0_wp
968                vanthoff                   = 2.31_wp
969             ENDIF
970          ENDIF
971
972!
973!--       Pre-calculate frequently calculated fractions of pi and rho_l
974          pirho_l  = pi * rho_l / 6.0_wp
975          dpirho_l = 1.0_wp / pirho_l
976
977          CALL location_message( 'finished', .TRUE. )
978
979       ELSE
980
981          CALL location_message( 'skipped', .TRUE. )
982
983       ENDIF
984
985    END SUBROUTINE bcm_init
986
987
988!------------------------------------------------------------------------------!
989! Description:
990! ------------
991!> Swapping of timelevels
992!------------------------------------------------------------------------------!
993    SUBROUTINE bcm_swap_timelevel ( mod_count )
994
995       IMPLICIT NONE
996
997       INTEGER, INTENT(IN) :: mod_count
998
999       IF ( bulk_cloud_model )  THEN
1000
1001          SELECT CASE ( mod_count )
1002
1003             CASE ( 0 )
1004
1005                IF ( microphysics_morrison )  THEN
1006                   qc => qc_1;    qc_p => qc_2
1007                   nc => nc_1;    nc_p => nc_2
1008                ENDIF
1009                IF ( microphysics_seifert )  THEN
1010                   qr => qr_1;    qr_p => qr_2
1011                   nr => nr_1;    nr_p => nr_2
1012                ENDIF
1013
1014             CASE ( 1 )
1015
1016                IF ( microphysics_morrison )  THEN
1017                   qc => qc_2;    qc_p => qc_1
1018                   nc => nc_2;    nc_p => nc_1
1019                ENDIF
1020                IF ( microphysics_seifert )  THEN
1021                   qr => qr_2;    qr_p => qr_1
1022                   nr => nr_2;    nr_p => nr_1
1023                ENDIF
1024
1025          END SELECT
1026
1027       ENDIF
1028
1029    END SUBROUTINE bcm_swap_timelevel
1030
1031
1032!------------------------------------------------------------------------------!
1033! Description:
1034! ------------
1035!> Header output for bulk cloud module
1036!------------------------------------------------------------------------------!
1037    SUBROUTINE bcm_header ( io )
1038
1039
1040       IMPLICIT NONE
1041
1042       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
1043
1044!
1045!--    Write bulk cloud module header
1046       WRITE ( io, 1 )
1047
1048       WRITE ( io, 2 )
1049       WRITE ( io, 3 )
1050
1051       IF ( microphysics_kessler )  THEN
1052          WRITE ( io, 4 ) 'Kessler-Scheme'
1053       ENDIF
1054
1055       IF ( microphysics_seifert )  THEN
1056          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
1057          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
1058          IF ( collision_turbulence )  WRITE ( io, 6 )
1059          IF ( ventilation_effect )  WRITE ( io, 7 )
1060          IF ( limiter_sedimentation )  WRITE ( io, 8 )
1061       ENDIF
1062
1063       WRITE ( io, 20 )
1064       WRITE ( io, 21 ) surface_pressure
1065       WRITE ( io, 22 ) r_d
1066       WRITE ( io, 23 ) rho_surface
1067       WRITE ( io, 24 ) c_p
1068       WRITE ( io, 25 ) l_v
1069
1070       IF ( microphysics_seifert )  THEN
1071          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
1072          WRITE ( io, 27 ) c_sedimentation
1073       ENDIF
1074
1075
10761   FORMAT ( //' Bulk cloud module information:'/ &
1077               ' ------------------------------------------'/ )
10782   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
1079             '    total water content is used.' )
10803   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
10814   FORMAT ( '--> Precipitation parameterization via ', A )
1082
10835   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
10846   FORMAT ( '--> Turbulence effects on precipitation process' )
10857   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
10868   FORMAT ( '--> Slope limiter used for sedimentation process' )
1087
108820  FORMAT ( '--> Essential parameters:' )
108921  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
109022  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
109123  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
109224  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
109325  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
109426  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
109527  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
1096
1097
1098    END SUBROUTINE bcm_header
1099
1100
1101!------------------------------------------------------------------------------!
1102!
1103! Description:
1104! ------------
1105!> Subroutine for averaging 3D data
1106!------------------------------------------------------------------------------!
1107    SUBROUTINE bcm_3d_data_averaging( mode, variable )
1108
1109       USE control_parameters,                                                 &
1110           ONLY:  average_count_3d
1111
1112       IMPLICIT NONE
1113
1114       CHARACTER (LEN=*) ::  mode     !<
1115       CHARACTER (LEN=*) ::  variable !<
1116
1117       INTEGER(iwp) ::  i       !< local index
1118       INTEGER(iwp) ::  j       !< local index
1119       INTEGER(iwp) ::  k       !< local index
1120
1121       IF ( mode == 'allocate' )  THEN
1122
1123          SELECT CASE ( TRIM( variable ) )
1124
1125             CASE ( 'nc' )
1126                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1127                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1128                ENDIF
1129                nc_av = 0.0_wp
1130
1131             CASE ( 'nr' )
1132                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1133                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1134                ENDIF
1135                nr_av = 0.0_wp
1136
1137             CASE ( 'prr' )
1138                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1139                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1140                ENDIF
1141                prr_av = 0.0_wp
1142
1143             CASE ( 'qc' )
1144                IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1145                   ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1146                ENDIF
1147                qc_av = 0.0_wp
1148
1149             CASE ( 'ql' )
1150                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1151                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1152                ENDIF
1153                ql_av = 0.0_wp
1154
1155             CASE ( 'qr' )
1156                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1157                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1158                ENDIF
1159                qr_av = 0.0_wp
1160
1161             CASE DEFAULT
1162                CONTINUE
1163
1164          END SELECT
1165
1166       ELSEIF ( mode == 'sum' )  THEN
1167
1168          SELECT CASE ( TRIM( variable ) )
1169
1170             CASE ( 'nc' )
1171                IF ( ALLOCATED( nc_av ) ) THEN
1172                   DO  i = nxlg, nxrg
1173                      DO  j = nysg, nyng
1174                         DO  k = nzb, nzt+1
1175                            nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
1176                         ENDDO
1177                      ENDDO
1178                   ENDDO
1179                ENDIF
1180
1181             CASE ( 'nr' )
1182                IF ( ALLOCATED( nr_av ) ) THEN
1183                   DO  i = nxlg, nxrg
1184                      DO  j = nysg, nyng
1185                         DO  k = nzb, nzt+1
1186                            nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
1187                         ENDDO
1188                      ENDDO
1189                   ENDDO
1190                ENDIF
1191
1192             CASE ( 'prr' )
1193                IF ( ALLOCATED( prr_av ) ) THEN
1194                   DO  i = nxlg, nxrg
1195                      DO  j = nysg, nyng
1196                         DO  k = nzb, nzt+1
1197                            prr_av(k,j,i) = prr_av(k,j,i) + prr(k,j,i)
1198                         ENDDO
1199                      ENDDO
1200                   ENDDO
1201                ENDIF
1202
1203             CASE ( 'qc' )
1204                IF ( ALLOCATED( qc_av ) ) THEN
1205                   DO  i = nxlg, nxrg
1206                      DO  j = nysg, nyng
1207                         DO  k = nzb, nzt+1
1208                            qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i)
1209                         ENDDO
1210                      ENDDO
1211                   ENDDO
1212                ENDIF
1213
1214             CASE ( 'ql' )
1215                IF ( ALLOCATED( ql_av ) ) THEN
1216                   DO  i = nxlg, nxrg
1217                      DO  j = nysg, nyng
1218                         DO  k = nzb, nzt+1
1219                            ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
1220                         ENDDO
1221                      ENDDO
1222                   ENDDO
1223                ENDIF
1224
1225             CASE ( 'qr' )
1226                IF ( ALLOCATED( qr_av ) ) THEN
1227                   DO  i = nxlg, nxrg
1228                      DO  j = nysg, nyng
1229                         DO  k = nzb, nzt+1
1230                            qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
1231                         ENDDO
1232                      ENDDO
1233                   ENDDO
1234                ENDIF
1235
1236             CASE DEFAULT
1237                CONTINUE
1238
1239          END SELECT
1240
1241       ELSEIF ( mode == 'average' )  THEN
1242
1243          SELECT CASE ( TRIM( variable ) )
1244
1245             CASE ( 'nc' )
1246                IF ( ALLOCATED( nc_av ) ) THEN
1247                   DO  i = nxlg, nxrg
1248                      DO  j = nysg, nyng
1249                         DO  k = nzb, nzt+1
1250                            nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1251                         ENDDO
1252                      ENDDO
1253                   ENDDO
1254                ENDIF
1255
1256             CASE ( 'nr' )
1257                IF ( ALLOCATED( nr_av ) ) THEN
1258                   DO  i = nxlg, nxrg
1259                      DO  j = nysg, nyng
1260                         DO  k = nzb, nzt+1
1261                            nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1262                         ENDDO
1263                      ENDDO
1264                   ENDDO
1265                ENDIF
1266
1267             CASE ( 'prr' )
1268                IF ( ALLOCATED( prr_av ) ) THEN
1269                   DO  i = nxlg, nxrg
1270                      DO  j = nysg, nyng
1271                         DO  k = nzb, nzt+1
1272                            prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1273                         ENDDO
1274                      ENDDO
1275                   ENDDO
1276                ENDIF
1277
1278             CASE ( 'qc' )
1279                IF ( ALLOCATED( qc_av ) ) THEN
1280                   DO  i = nxlg, nxrg
1281                      DO  j = nysg, nyng
1282                         DO  k = nzb, nzt+1
1283                            qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1284                         ENDDO
1285                      ENDDO
1286                   ENDDO
1287                ENDIF
1288
1289             CASE ( 'ql' )
1290                IF ( ALLOCATED( ql_av ) ) THEN
1291                   DO  i = nxlg, nxrg
1292                      DO  j = nysg, nyng
1293                         DO  k = nzb, nzt+1
1294                            ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1295                         ENDDO
1296                      ENDDO
1297                   ENDDO
1298                ENDIF
1299
1300             CASE ( 'qr' )
1301                IF ( ALLOCATED( qr_av ) ) THEN
1302                   DO  i = nxlg, nxrg
1303                      DO  j = nysg, nyng
1304                         DO  k = nzb, nzt+1
1305                            qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1306                         ENDDO
1307                      ENDDO
1308                   ENDDO
1309                ENDIF
1310
1311             CASE DEFAULT
1312                CONTINUE
1313
1314          END SELECT
1315
1316       ENDIF
1317
1318    END SUBROUTINE bcm_3d_data_averaging
1319
1320
1321!------------------------------------------------------------------------------!
1322! Description:
1323! ------------
1324!> Define 2D output variables.
1325!------------------------------------------------------------------------------!
1326 SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
1327                                two_d, nzb_do, nzt_do )
1328
1329
1330    IMPLICIT NONE
1331
1332    CHARACTER (LEN=*), INTENT(INOUT) ::  grid       !< name of vertical grid
1333    CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
1334    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1335
1336    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1337    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< vertical output index (bottom)
1338    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical output index (top)
1339
1340    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1341
1342    INTEGER(iwp) ::  i         !< loop index along x-direction
1343    INTEGER(iwp) ::  j         !< loop index along y-direction
1344    INTEGER(iwp) ::  k         !< loop index along z-direction
1345
1346    LOGICAL, INTENT(INOUT) ::  found   !< flag if output variable is found
1347    LOGICAL, INTENT(INOUT) ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1348    LOGICAL ::  resorted  !< flag if output is already resorted
1349
1350    REAL(wp), PARAMETER ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1351
1352    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf !< local
1353       !< array to which output data is resorted to
1354
1355    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1356
1357    found = .TRUE.
1358    resorted = .FALSE.
1359!
1360!-- Set masking flag for topography for not resorted arrays
1361    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1362
1363    SELECT CASE ( TRIM( variable ) )
1364
1365       CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
1366          IF ( av == 0 )  THEN
1367             to_be_resorted => nc
1368          ELSE
1369             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1370                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1371                nc_av = REAL( fill_value, KIND = wp )
1372             ENDIF
1373             to_be_resorted => nc_av
1374          ENDIF
1375          IF ( mode == 'xy' ) grid = 'zu'
1376
1377       CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
1378          IF ( av == 0 )  THEN
1379             to_be_resorted => nr
1380          ELSE
1381             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1382                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1383                nr_av = REAL( fill_value, KIND = wp )
1384             ENDIF
1385             to_be_resorted => nr_av
1386          ENDIF
1387          IF ( mode == 'xy' ) grid = 'zu'
1388
1389       CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
1390!                CALL exchange_horiz_2d( precipitation_amount )
1391          DO  i = nxl, nxr
1392             DO  j = nys, nyn
1393                local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
1394             ENDDO
1395          ENDDO
1396          precipitation_amount = 0.0_wp   ! reset for next integ. interval
1397          resorted = .TRUE.
1398          two_d = .TRUE.
1399          IF ( mode == 'xy' ) grid = 'zu1'
1400
1401       CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
1402          IF ( av == 0 )  THEN
1403!                   CALL exchange_horiz( prr, nbgp )
1404             DO  i = nxl, nxr
1405                DO  j = nys, nyn
1406                   DO  k = nzb, nzt+1
1407                      local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
1408                   ENDDO
1409                ENDDO
1410             ENDDO
1411          ELSE
1412             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1413                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1414                prr_av = REAL( fill_value, KIND = wp )
1415             ENDIF
1416!                   CALL exchange_horiz( prr_av, nbgp )
1417             DO  i = nxl, nxr
1418                DO  j = nys, nyn
1419                   DO  k = nzb, nzt+1
1420                      local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
1421                   ENDDO
1422                ENDDO
1423             ENDDO
1424          ENDIF
1425          resorted = .TRUE.
1426          IF ( mode == 'xy' ) grid = 'zu'
1427
1428       CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
1429          IF ( av == 0 )  THEN
1430             to_be_resorted => qc
1431          ELSE
1432             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1433                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1434                qc_av = REAL( fill_value, KIND = wp )
1435             ENDIF
1436             to_be_resorted => qc_av
1437          ENDIF
1438          IF ( mode == 'xy' ) grid = 'zu'
1439
1440       CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
1441          IF ( av == 0 )  THEN
1442             to_be_resorted => ql
1443          ELSE
1444             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1445                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1446                ql_av = REAL( fill_value, KIND = wp )
1447             ENDIF
1448             to_be_resorted => ql_av
1449          ENDIF
1450          IF ( mode == 'xy' ) grid = 'zu'
1451
1452       CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
1453          IF ( av == 0 )  THEN
1454             to_be_resorted => qr
1455          ELSE
1456             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1457                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1458                qr_av = REAL( fill_value, KIND = wp )
1459             ENDIF
1460             to_be_resorted => qr_av
1461          ENDIF
1462          IF ( mode == 'xy' ) grid = 'zu'
1463
1464       CASE DEFAULT
1465          found = .FALSE.
1466          grid  = 'none'
1467
1468    END SELECT
1469
1470    IF ( found .AND. .NOT. resorted )  THEN
1471       DO  i = nxl, nxr
1472          DO  j = nys, nyn
1473             DO  k = nzb_do, nzt_do
1474                local_pf(i,j,k) = MERGE( &
1475                                     to_be_resorted(k,j,i),                    &
1476                                     REAL( fill_value, KIND = wp ),            &
1477                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1478                                  )
1479             ENDDO
1480          ENDDO
1481       ENDDO
1482    ENDIF
1483
1484 END SUBROUTINE bcm_data_output_2d
1485
1486
1487!------------------------------------------------------------------------------!
1488! Description:
1489! ------------
1490!> Define 3D output variables.
1491!------------------------------------------------------------------------------!
1492 SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1493
1494
1495    IMPLICIT NONE
1496
1497    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1498
1499    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1500    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< lower limit of the data output (usually 0)
1501    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
1502
1503    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1504
1505    INTEGER(iwp) ::  i         !< loop index along x-direction
1506    INTEGER(iwp) ::  j         !< loop index along y-direction
1507    INTEGER(iwp) ::  k         !< loop index along z-direction
1508
1509    LOGICAL, INTENT(INOUT) ::  found     !< flag if output variable is found
1510    LOGICAL ::  resorted  !< flag if output is already resorted
1511
1512    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1513
1514    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf   !< local
1515       !< array to which output data is resorted to
1516
1517    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1518
1519    found = .TRUE.
1520    resorted = .FALSE.
1521!
1522!-- Set masking flag for topography for not resorted arrays
1523    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1524
1525    SELECT CASE ( TRIM( variable ) )
1526
1527       CASE ( 'nc' )
1528          IF ( av == 0 )  THEN
1529             to_be_resorted => nc
1530          ELSE
1531             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1532                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1533                nc_av = REAL( fill_value, KIND = wp )
1534             ENDIF
1535             to_be_resorted => nc_av
1536          ENDIF
1537
1538       CASE ( 'nr' )
1539          IF ( av == 0 )  THEN
1540             to_be_resorted => nr
1541          ELSE
1542             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1543                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1544                nr_av = REAL( fill_value, KIND = wp )
1545             ENDIF
1546             to_be_resorted => nr_av
1547          ENDIF
1548
1549       CASE ( 'prr' )
1550          IF ( av == 0 )  THEN
1551             DO  i = nxl, nxr
1552                DO  j = nys, nyn
1553                   DO  k = nzb_do, nzt_do
1554                      local_pf(i,j,k) = prr(k,j,i)
1555                   ENDDO
1556                ENDDO
1557             ENDDO
1558          ELSE
1559             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1560                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1561                prr_av = REAL( fill_value, KIND = wp )
1562             ENDIF
1563             DO  i = nxl, nxr
1564                DO  j = nys, nyn
1565                   DO  k = nzb_do, nzt_do
1566                      local_pf(i,j,k) = prr_av(k,j,i)
1567                   ENDDO
1568                ENDDO
1569             ENDDO
1570          ENDIF
1571          resorted = .TRUE.
1572
1573       CASE ( 'qc' )
1574          IF ( av == 0 )  THEN
1575             to_be_resorted => qc
1576          ELSE
1577             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1578                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1579                qc_av = REAL( fill_value, KIND = wp )
1580             ENDIF
1581             to_be_resorted => qc_av
1582          ENDIF
1583
1584       CASE ( 'ql' )
1585          IF ( av == 0 )  THEN
1586             to_be_resorted => ql
1587          ELSE
1588             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1589                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1590                ql_av = REAL( fill_value, KIND = wp )
1591             ENDIF
1592             to_be_resorted => ql_av
1593          ENDIF
1594
1595       CASE ( 'qr' )
1596          IF ( av == 0 )  THEN
1597             to_be_resorted => qr
1598          ELSE
1599             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1600                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1601                qr_av = REAL( fill_value, KIND = wp )
1602             ENDIF
1603             to_be_resorted => qr_av
1604          ENDIF
1605
1606       CASE DEFAULT
1607          found = .FALSE.
1608
1609    END SELECT
1610
1611
1612    IF ( found .AND. .NOT. resorted )  THEN
1613       DO  i = nxl, nxr
1614          DO  j = nys, nyn
1615             DO  k = nzb_do, nzt_do
1616                local_pf(i,j,k) = MERGE(                                       &
1617                                     to_be_resorted(k,j,i),                    &
1618                                     REAL( fill_value, KIND = wp ),            &
1619                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1620                                  )
1621             ENDDO
1622          ENDDO
1623       ENDDO
1624    ENDIF
1625
1626 END SUBROUTINE bcm_data_output_3d
1627
1628
1629!------------------------------------------------------------------------------!
1630! Description:
1631! ------------
1632!> This routine reads the respective restart data for the bulk cloud module.
1633!------------------------------------------------------------------------------!
1634    SUBROUTINE bcm_rrd_global( found )
1635
1636
1637       USE control_parameters,                                                 &
1638           ONLY: length, restart_string
1639
1640
1641       IMPLICIT NONE
1642
1643       LOGICAL, INTENT(OUT)  ::  found
1644
1645
1646       found = .TRUE.
1647
1648       SELECT CASE ( restart_string(1:length) )
1649
1650          CASE ( 'c_sedimentation' )
1651             READ ( 13 )  c_sedimentation
1652
1653          CASE ( 'bulk_cloud_model' )
1654             READ ( 13 )  bulk_cloud_model
1655
1656          CASE ( 'cloud_scheme' )
1657             READ ( 13 )  cloud_scheme
1658
1659          CASE ( 'cloud_water_sedimentation' )
1660             READ ( 13 )  cloud_water_sedimentation
1661
1662          CASE ( 'collision_turbulence' )
1663             READ ( 13 )  collision_turbulence
1664
1665          CASE ( 'limiter_sedimentation' )
1666             READ ( 13 )  limiter_sedimentation
1667
1668          CASE ( 'nc_const' )
1669             READ ( 13 )  nc_const
1670
1671          CASE ( 'precipitation' )
1672             READ ( 13 ) precipitation
1673
1674          CASE ( 'ventilation_effect' )
1675             READ ( 13 )  ventilation_effect
1676
1677          CASE ( 'na_init' )
1678             READ ( 13 )  na_init
1679
1680          CASE ( 'dry_aerosol_radius' )
1681             READ ( 13 )  dry_aerosol_radius
1682
1683          CASE ( 'sigma_bulk' )
1684             READ ( 13 )  sigma_bulk
1685
1686          CASE ( 'aerosol_bulk' )
1687             READ ( 13 )  aerosol_bulk
1688
1689          CASE ( 'curvature_solution_effects_bulk' )
1690             READ ( 13 )  curvature_solution_effects_bulk
1691
1692
1693!          CASE ( 'global_paramter' )
1694!             READ ( 13 )  global_parameter
1695!          CASE ( 'global_array' )
1696!             IF ( .NOT. ALLOCATED( global_array ) )  ALLOCATE( global_array(1:10) )
1697!             READ ( 13 )  global_array
1698
1699          CASE DEFAULT
1700
1701             found = .FALSE.
1702
1703       END SELECT
1704
1705
1706    END SUBROUTINE bcm_rrd_global
1707
1708
1709!------------------------------------------------------------------------------!
1710! Description:
1711! ------------
1712!> This routine reads the respective restart data for the bulk cloud module.
1713!------------------------------------------------------------------------------!
1714    SUBROUTINE bcm_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,          &
1715                              nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
1716                              nysc, nys_on_file, tmp_2d, tmp_3d, found )
1717
1718
1719       USE control_parameters
1720
1721       USE indices
1722
1723       USE pegrid
1724
1725
1726       IMPLICIT NONE
1727
1728       INTEGER(iwp) ::  k               !<
1729       INTEGER(iwp) ::  nxlc            !<
1730       INTEGER(iwp) ::  nxlf            !<
1731       INTEGER(iwp) ::  nxl_on_file     !<
1732       INTEGER(iwp) ::  nxrc            !<
1733       INTEGER(iwp) ::  nxrf            !<
1734       INTEGER(iwp) ::  nxr_on_file     !<
1735       INTEGER(iwp) ::  nync            !<
1736       INTEGER(iwp) ::  nynf            !<
1737       INTEGER(iwp) ::  nyn_on_file     !<
1738       INTEGER(iwp) ::  nysc            !<
1739       INTEGER(iwp) ::  nysf            !<
1740       INTEGER(iwp) ::  nys_on_file     !<
1741
1742       LOGICAL, INTENT(OUT)  ::  found
1743
1744       REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1745       REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1746
1747!
1748!-- Here the reading of user-defined restart data follows:
1749!-- Sample for user-defined output
1750
1751
1752       found = .TRUE.
1753
1754       SELECT CASE ( restart_string(1:length) )
1755
1756          CASE ( 'nc' )
1757             IF ( k == 1 )  READ ( 13 )  tmp_3d
1758             nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1759                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1760
1761          CASE ( 'nc_av' )
1762             IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1763                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1764             ENDIF
1765             IF ( k == 1 )  READ ( 13 )  tmp_3d
1766             nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1767                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1768
1769          CASE ( 'nr' )
1770             IF ( k == 1 )  READ ( 13 )  tmp_3d
1771             nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1772                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1773
1774          CASE ( 'nr_av' )
1775             IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1776                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1777             ENDIF
1778             IF ( k == 1 )  READ ( 13 )  tmp_3d
1779             nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1780                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1781
1782          CASE ( 'precipitation_amount' )
1783             IF ( k == 1 )  READ ( 13 )  tmp_2d
1784             precipitation_amount(nysc-nbgp:nync+nbgp,                   &
1785                                  nxlc-nbgp:nxrc+nbgp)  =                &
1786                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1787
1788          CASE ( 'prr' )
1789             IF ( .NOT. ALLOCATED( prr ) )  THEN
1790                ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1791             ENDIF
1792             IF ( k == 1 )  READ ( 13 )  tmp_3d
1793             prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
1794                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1795
1796          CASE ( 'prr_av' )
1797             IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1798                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1799             ENDIF
1800             IF ( k == 1 )  READ ( 13 )  tmp_3d
1801             prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
1802                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1803
1804          CASE ( 'qc' )
1805             IF ( k == 1 )  READ ( 13 )  tmp_3d
1806             qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1807                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1808
1809          CASE ( 'qc_av' )
1810             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1811                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1812             ENDIF
1813             IF ( k == 1 )  READ ( 13 )  tmp_3d
1814             qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1815                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1816
1817          CASE ( 'ql' )
1818             IF ( k == 1 )  READ ( 13 )  tmp_3d
1819             ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1820                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1821
1822          CASE ( 'ql_av' )
1823             IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1824                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1825             ENDIF
1826             IF ( k == 1 )  READ ( 13 )  tmp_3d
1827             ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1828                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1829
1830          CASE ( 'qr' )
1831             IF ( k == 1 )  READ ( 13 )  tmp_3d
1832             qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1833                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1834
1835          CASE ( 'qr_av' )
1836             IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1837                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1838             ENDIF
1839             IF ( k == 1 )  READ ( 13 )  tmp_3d
1840             qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1841                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1842!
1843          CASE DEFAULT
1844
1845             found = .FALSE.
1846
1847          END SELECT
1848
1849
1850    END SUBROUTINE bcm_rrd_local
1851
1852
1853!------------------------------------------------------------------------------!
1854! Description:
1855! ------------
1856!> This routine writes the respective restart data for the bulk cloud module.
1857!------------------------------------------------------------------------------!
1858    SUBROUTINE bcm_wrd_global
1859
1860
1861       IMPLICIT NONE
1862
1863       CALL wrd_write_string( 'c_sedimentation' )
1864       WRITE ( 14 )  c_sedimentation
1865
1866       CALL wrd_write_string( 'bulk_cloud_model' )
1867       WRITE ( 14 )  bulk_cloud_model
1868
1869       CALL wrd_write_string( 'cloud_scheme' )
1870       WRITE ( 14 )  cloud_scheme
1871
1872       CALL wrd_write_string( 'cloud_water_sedimentation' )
1873       WRITE ( 14 )  cloud_water_sedimentation
1874
1875       CALL wrd_write_string( 'collision_turbulence' )
1876       WRITE ( 14 )  collision_turbulence
1877
1878       CALL wrd_write_string( 'limiter_sedimentation' )
1879       WRITE ( 14 )  limiter_sedimentation
1880
1881       CALL wrd_write_string( 'nc_const' )
1882       WRITE ( 14 )  nc_const
1883
1884       CALL wrd_write_string( 'precipitation' )
1885       WRITE ( 14 )  precipitation
1886
1887       CALL wrd_write_string( 'ventilation_effect' )
1888       WRITE ( 14 )  ventilation_effect
1889
1890       CALL wrd_write_string( 'na_init' )
1891       WRITE ( 14 )  na_init
1892
1893       CALL wrd_write_string( 'dry_aerosol_radius' )
1894       WRITE ( 14 )  dry_aerosol_radius
1895
1896       CALL wrd_write_string( 'sigma_bulk' )
1897       WRITE ( 14 )  sigma_bulk
1898
1899       CALL wrd_write_string( 'aerosol_bulk' )
1900       WRITE ( 14 )  aerosol_bulk
1901
1902       CALL wrd_write_string( 'curvature_solution_effects_bulk' )
1903       WRITE ( 14 )  curvature_solution_effects_bulk
1904
1905
1906! needs preceeding allocation if array
1907!       CALL wrd_write_string( 'global_parameter' )
1908!       WRITE ( 14 )  global_parameter
1909
1910!       IF ( ALLOCATED( inflow_damping_factor ) )  THEN
1911!          CALL wrd_write_string( 'inflow_damping_factor' )
1912!          WRITE ( 14 )  inflow_damping_factor
1913!       ENDIF
1914
1915
1916    END SUBROUTINE bcm_wrd_global
1917
1918
1919!------------------------------------------------------------------------------!
1920! Description:
1921! ------------
1922!> This routine writes the respective restart data for the bulk cloud module.
1923!------------------------------------------------------------------------------!
1924    SUBROUTINE bcm_wrd_local
1925
1926
1927       IMPLICIT NONE
1928
1929       IF ( ALLOCATED( prr ) )  THEN
1930          CALL wrd_write_string( 'prr' )
1931          WRITE ( 14 )  prr
1932       ENDIF
1933
1934       IF ( ALLOCATED( prr_av ) )  THEN
1935          CALL wrd_write_string( 'prr_av' )
1936          WRITE ( 14 )  prr_av
1937       ENDIF
1938
1939       IF ( ALLOCATED( precipitation_amount ) )  THEN
1940          CALL wrd_write_string( 'precipitation_amount' )
1941          WRITE ( 14 )  precipitation_amount
1942       ENDIF
1943
1944       CALL wrd_write_string( 'ql' )
1945       WRITE ( 14 )  ql
1946
1947       IF ( ALLOCATED( ql_av ) )  THEN
1948          CALL wrd_write_string( 'ql_av' )
1949          WRITE ( 14 )  ql_av
1950       ENDIF
1951
1952       CALL wrd_write_string( 'qc' )
1953       WRITE ( 14 )  qc
1954
1955       IF ( ALLOCATED( qc_av ) )  THEN
1956          CALL wrd_write_string( 'qc_av' )
1957          WRITE ( 14 )  qc_av
1958       ENDIF
1959
1960       IF ( microphysics_morrison )  THEN
1961
1962          CALL wrd_write_string( 'nc' )
1963          WRITE ( 14 )  nc
1964
1965          IF ( ALLOCATED( nc_av ) )  THEN
1966             CALL wrd_write_string( 'nc_av' )
1967             WRITE ( 14 )  nc_av
1968          ENDIF
1969
1970       ENDIF
1971
1972       IF ( microphysics_seifert )  THEN
1973
1974          CALL wrd_write_string( 'nr' )
1975          WRITE ( 14 )  nr
1976
1977          IF ( ALLOCATED( nr_av ) )  THEN
1978             CALL wrd_write_string( 'nr_av' )
1979             WRITE ( 14 )  nr_av
1980          ENDIF
1981
1982          CALL wrd_write_string( 'qr' )
1983          WRITE ( 14 )  qr
1984
1985          IF ( ALLOCATED( qr_av ) )  THEN
1986             CALL wrd_write_string( 'qr_av' )
1987             WRITE ( 14 )  qr_av
1988          ENDIF
1989
1990       ENDIF
1991
1992
1993    END SUBROUTINE bcm_wrd_local
1994
1995
1996!------------------------------------------------------------------------------!
1997! Description:
1998! ------------
1999!> Control of microphysics for all grid points
2000!------------------------------------------------------------------------------!
2001    SUBROUTINE bcm_actions
2002
2003       IMPLICIT NONE
2004
2005       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
2006!
2007!--       Calculate vertical profile of the hydrostatic pressure (hyp)
2008          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
2009          d_exner = exner_function_invers(hyp)
2010          exner = 1.0_wp / exner_function_invers(hyp)
2011          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
2012!
2013!--       Compute reference density
2014          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
2015       ENDIF
2016
2017!
2018!--    Compute length of time step
2019       IF ( call_microphysics_at_all_substeps )  THEN
2020          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
2021       ELSE
2022          dt_micro = dt_3d
2023       ENDIF
2024
2025!
2026!--    Reset precipitation rate
2027       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
2028
2029!
2030!--    Compute cloud physics
2031       IF ( microphysics_kessler )  THEN
2032
2033          CALL autoconversion_kessler
2034          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2035
2036       ELSEIF ( microphysics_seifert )  THEN
2037
2038          CALL adjust_cloud
2039          IF ( microphysics_morrison )  CALL activation
2040          IF ( microphysics_morrison )  CALL condensation
2041          CALL autoconversion
2042          CALL accretion
2043          CALL selfcollection_breakup
2044          CALL evaporation_rain
2045          CALL sedimentation_rain
2046          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2047
2048       ENDIF
2049
2050       CALL calc_precipitation_amount
2051
2052    END SUBROUTINE bcm_actions
2053
2054!------------------------------------------------------------------------------!
2055! Description:
2056! ------------
2057!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
2058!> evaporation of rain drops due to too small or too big weights
2059!> of rain drops (Stevens and Seifert, 2008).
2060!------------------------------------------------------------------------------!
2061    SUBROUTINE adjust_cloud
2062
2063       IMPLICIT NONE
2064
2065       INTEGER(iwp) ::  i                 !<
2066       INTEGER(iwp) ::  j                 !<
2067       INTEGER(iwp) ::  k                 !<
2068
2069       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2070
2071       CALL cpu_log( log_point_s(50), 'adjust_cloud', 'start' )
2072
2073       DO  i = nxlg, nxrg
2074          DO  j = nysg, nyng
2075             DO  k = nzb+1, nzt
2076!
2077!--             Predetermine flag to mask topography
2078                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2079
2080                IF ( qr(k,j,i) <= eps_sb )  THEN
2081                   qr(k,j,i) = 0.0_wp
2082                   nr(k,j,i) = 0.0_wp
2083                ELSE
2084                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
2085                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
2086                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
2087                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
2088                   ENDIF
2089                ENDIF
2090
2091                IF ( microphysics_morrison ) THEN
2092                   IF ( qc(k,j,i) <= eps_sb )  THEN
2093                      qc(k,j,i) = 0.0_wp
2094                      nc(k,j,i) = 0.0_wp
2095                   ELSE
2096                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2097                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * flag
2098                      ENDIF
2099                   ENDIF
2100                ENDIF
2101
2102             ENDDO
2103          ENDDO
2104       ENDDO
2105
2106       CALL cpu_log( log_point_s(50), 'adjust_cloud', 'stop' )
2107
2108    END SUBROUTINE adjust_cloud
2109
2110!------------------------------------------------------------------------------!
2111! Description:
2112! ------------
2113!> Calculate number of activated condensation nucleii after simple activation
2114!> scheme of Twomey, 1959.
2115!------------------------------------------------------------------------------!
2116    SUBROUTINE activation
2117
2118       IMPLICIT NONE
2119
2120       INTEGER(iwp) ::  i                 !<
2121       INTEGER(iwp) ::  j                 !<
2122       INTEGER(iwp) ::  k                 !<
2123
2124       REAL(wp)     ::  activ             !<
2125       REAL(wp)     ::  afactor           !<
2126       REAL(wp)     ::  beta_act          !<
2127       REAL(wp)     ::  bfactor           !<
2128       REAL(wp)     ::  k_act             !<
2129       REAL(wp)     ::  n_act             !<
2130       REAL(wp)     ::  n_ccn             !<
2131       REAL(wp)     ::  s_0               !<
2132       REAL(wp)     ::  sat_max           !<
2133       REAL(wp)     ::  sigma             !<
2134       REAL(wp)     ::  sigma_act         !<
2135
2136       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2137
2138       CALL cpu_log( log_point_s(65), 'activation', 'start' )
2139
2140       DO  i = nxlg, nxrg
2141          DO  j = nysg, nyng
2142             DO  k = nzb+1, nzt
2143!
2144!--             Predetermine flag to mask topography
2145                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2146
2147!
2148!--             Call calculation of supersaturation located in subroutine
2149                CALL supersaturation ( i, j, k )
2150!
2151!--             Prescribe parameters for activation
2152!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
2153                k_act  = 0.7_wp
2154                activ  = 0.0_wp
2155
2156
2157                IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
2158!
2159!--                Compute the number of activated Aerosols
2160!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
2161                   n_act     = na_init * sat**k_act
2162!
2163!--                Compute the number of cloud droplets
2164!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2165!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
2166
2167!
2168!--                Compute activation rate after Khairoutdinov and Kogan
2169!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2170                   sat_max = 1.0_wp / 100.0_wp
2171                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
2172                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
2173                       dt_micro
2174                ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
2175!
2176!--                Curvature effect (afactor) with surface tension
2177!--                parameterization by Straka (2009)
2178                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
2179                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
2180!
2181!--                Solute effect (bfactor)
2182                   bfactor = vanthoff * molecular_weight_of_water *            &
2183                       rho_s / ( molecular_weight_of_solute * rho_l )
2184
2185!
2186!--                Prescribe power index that describes the soluble fraction
2187!--                of an aerosol particle (beta)
2188!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2189                   beta_act  = 0.5_wp
2190                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2191!
2192!--                Calculate mean geometric supersaturation (s_0) with
2193!--                parameterization by Khvorostyanov and Curry (2006)
2194                   s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
2195                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
2196
2197!
2198!--                Calculate number of activated CCN as a function of
2199!--                supersaturation and taking Koehler theory into account
2200!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2201                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
2202                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2203                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
2204                ENDIF
2205
2206                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
2207
2208             ENDDO
2209          ENDDO
2210       ENDDO
2211
2212       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
2213
2214    END SUBROUTINE activation
2215
2216
2217!------------------------------------------------------------------------------!
2218! Description:
2219! ------------
2220!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2221!> Kogan, 2000).
2222!------------------------------------------------------------------------------!
2223    SUBROUTINE condensation
2224
2225       IMPLICIT NONE
2226
2227       INTEGER(iwp) ::  i                 !<
2228       INTEGER(iwp) ::  j                 !<
2229       INTEGER(iwp) ::  k                 !<
2230
2231       REAL(wp)     ::  cond              !<
2232       REAL(wp)     ::  cond_max          !<
2233       REAL(wp)     ::  dc                !<
2234       REAL(wp)     ::  evap              !<
2235       REAL(wp)     ::  g_fac             !<
2236       REAL(wp)     ::  nc_0              !<
2237       REAL(wp)     ::  temp              !<
2238       REAL(wp)     ::  xc                !<
2239
2240       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2241
2242       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
2243
2244       DO  i = nxlg, nxrg
2245          DO  j = nysg, nyng
2246             DO  k = nzb+1, nzt
2247!
2248!--             Predetermine flag to mask topography
2249                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2250!
2251!--             Call calculation of supersaturation
2252                CALL supersaturation ( i, j, k )
2253!
2254!--             Actual temperature:
2255                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2256
2257                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
2258                                    l_v / ( thermal_conductivity_l * temp )    &
2259                                    + r_v * temp / ( diff_coeff_l * e_s )      &
2260                                    )
2261!
2262!--             Mean weight of cloud drops
2263                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
2264                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
2265!
2266!--             Weight averaged diameter of cloud drops:
2267                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2268!
2269!--             Integral diameter of cloud drops
2270                nc_0 = nc(k,j,i) * dc
2271!
2272!--             Condensation needs only to be calculated in supersaturated regions
2273                IF ( sat > 0.0_wp )  THEN
2274!
2275!--                Condensation rate of cloud water content
2276!--                after KK scheme.
2277!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2278                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2279                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
2280                   cond      = MIN( cond, cond_max / dt_micro )
2281
2282                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
2283                ELSEIF ( sat < 0.0_wp ) THEN
2284                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2285                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
2286
2287                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
2288                ENDIF
2289                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2290                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
2291                ENDIF
2292             ENDDO
2293          ENDDO
2294       ENDDO
2295
2296       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
2297
2298    END SUBROUTINE condensation
2299
2300
2301!------------------------------------------------------------------------------!
2302! Description:
2303! ------------
2304!> Autoconversion rate (Seifert and Beheng, 2006).
2305!------------------------------------------------------------------------------!
2306    SUBROUTINE autoconversion
2307
2308       IMPLICIT NONE
2309
2310       INTEGER(iwp) ::  i                 !<
2311       INTEGER(iwp) ::  j                 !<
2312       INTEGER(iwp) ::  k                 !<
2313
2314       REAL(wp)     ::  alpha_cc          !<
2315       REAL(wp)     ::  autocon           !<
2316       REAL(wp)     ::  dissipation       !<
2317       REAL(wp)     ::  flag              !< flag to mask topography grid points
2318       REAL(wp)     ::  k_au              !<
2319       REAL(wp)     ::  l_mix             !<
2320       REAL(wp)     ::  nc_auto           !<
2321       REAL(wp)     ::  nu_c              !<
2322       REAL(wp)     ::  phi_au            !<
2323       REAL(wp)     ::  r_cc              !<
2324       REAL(wp)     ::  rc                !<
2325       REAL(wp)     ::  re_lambda         !<
2326       REAL(wp)     ::  sigma_cc          !<
2327       REAL(wp)     ::  tau_cloud         !<
2328       REAL(wp)     ::  xc                !<
2329
2330       CALL cpu_log( log_point_s(47), 'autoconversion', 'start' )
2331
2332       DO  i = nxlg, nxrg
2333          DO  j = nysg, nyng
2334             DO  k = nzb+1, nzt
2335!
2336!--             Predetermine flag to mask topography
2337                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2338
2339                IF ( microphysics_morrison ) THEN
2340                   nc_auto = nc(k,j,i)
2341                ELSE
2342                   nc_auto = nc_const
2343                ENDIF
2344
2345                IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
2346
2347                   k_au = k_cc / ( 20.0_wp * x0 )
2348!
2349!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2350!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
2351                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
2352                                    qc(k,j,i) ), 0.0_wp )
2353!
2354!--                Universal function for autoconversion process
2355!--                (Seifert and Beheng, 2006):
2356                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
2357                            ( 1.0_wp - tau_cloud**0.68_wp )**3
2358!
2359!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
2360!--                (Use constant nu_c = 1.0_wp instead?)
2361                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
2362!
2363!--                Mean weight of cloud droplets:
2364                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
2365!
2366!--                Parameterized turbulence effects on autoconversion (Seifert,
2367!--                Nuijens and Stevens, 2010)
2368                   IF ( collision_turbulence )  THEN
2369!
2370!--                   Weight averaged radius of cloud droplets:
2371                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2372
2373                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2374                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2375                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2376!
2377!--                   Mixing length (neglecting distance to ground and
2378!--                   stratification)
2379                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2380!
2381!--                   Limit dissipation rate according to Seifert, Nuijens and
2382!--                   Stevens (2010)
2383                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
2384!
2385!--                   Compute Taylor-microscale Reynolds number:
2386                      re_lambda = 6.0_wp / 11.0_wp *                           &
2387                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
2388                                  SQRT( 15.0_wp / kin_vis_air ) *              &
2389                                  dissipation**( 1.0_wp / 6.0_wp )
2390!
2391!--                   The factor of 1.0E4 is needed to convert the dissipation
2392!--                   rate from m2 s-3 to cm2 s-3.
2393                      k_au = k_au * ( 1.0_wp +                                 &
2394                                      dissipation * 1.0E4_wp *                 &
2395                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
2396                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
2397                                                                      r_cc ) / &
2398                                                        sigma_cc )**2          &
2399                                                      ) + beta_cc              &
2400                                      )                                        &
2401                                    )
2402                   ENDIF
2403!
2404!--                Autoconversion rate (Seifert and Beheng, 2006):
2405                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
2406                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
2407                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
2408                             rho_surface
2409                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
2410
2411                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
2412                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
2413                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
2414                                                              * flag
2415                   IF ( microphysics_morrison ) THEN
2416                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
2417                                    autocon / x0 * hyrho(k) * dt_micro * flag )
2418                   ENDIF
2419
2420                ENDIF
2421
2422             ENDDO
2423          ENDDO
2424       ENDDO
2425
2426       CALL cpu_log( log_point_s(47), 'autoconversion', 'stop' )
2427
2428    END SUBROUTINE autoconversion
2429
2430
2431!------------------------------------------------------------------------------!
2432! Description:
2433! ------------
2434!> Autoconversion process (Kessler, 1969).
2435!------------------------------------------------------------------------------!
2436    SUBROUTINE autoconversion_kessler
2437
2438
2439       IMPLICIT NONE
2440
2441       INTEGER(iwp) ::  i      !<
2442       INTEGER(iwp) ::  j      !<
2443       INTEGER(iwp) ::  k      !<
2444       INTEGER(iwp) ::  k_wall !< topgraphy top index
2445
2446       REAL(wp)    ::  dqdt_precip !<
2447       REAL(wp)    ::  flag        !< flag to mask topography grid points
2448
2449       DO  i = nxlg, nxrg
2450          DO  j = nysg, nyng
2451!
2452!--          Determine vertical index of topography top
2453             k_wall = get_topography_top_index_ji( j, i, 's' )
2454             DO  k = nzb+1, nzt
2455!
2456!--             Predetermine flag to mask topography
2457                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2458
2459                IF ( qc(k,j,i) > ql_crit )  THEN
2460                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
2461                ELSE
2462                   dqdt_precip = 0.0_wp
2463                ENDIF
2464
2465                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
2466                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
2467                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp *     &
2468                                        d_exner(k)              * flag
2469
2470!
2471!--             Compute the rain rate (stored on surface grid point)
2472                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2473
2474             ENDDO
2475          ENDDO
2476       ENDDO
2477
2478   END SUBROUTINE autoconversion_kessler
2479
2480
2481!------------------------------------------------------------------------------!
2482! Description:
2483! ------------
2484!> Accretion rate (Seifert and Beheng, 2006).
2485!------------------------------------------------------------------------------!
2486    SUBROUTINE accretion
2487
2488       IMPLICIT NONE
2489
2490       INTEGER(iwp) ::  i                 !<
2491       INTEGER(iwp) ::  j                 !<
2492       INTEGER(iwp) ::  k                 !<
2493
2494       REAL(wp)     ::  accr              !<
2495       REAL(wp)     ::  flag              !< flag to mask topography grid points
2496       REAL(wp)     ::  k_cr              !<
2497       REAL(wp)     ::  nc_accr           !<
2498       REAL(wp)     ::  phi_ac            !<
2499       REAL(wp)     ::  tau_cloud         !<
2500       REAL(wp)     ::  xc                !<
2501
2502
2503       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
2504
2505       DO  i = nxlg, nxrg
2506          DO  j = nysg, nyng
2507             DO  k = nzb+1, nzt
2508!
2509!--             Predetermine flag to mask topography
2510                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2511
2512                IF ( microphysics_morrison ) THEN
2513                   nc_accr = nc(k,j,i)
2514                ELSE
2515                   nc_accr = nc_const
2516                ENDIF
2517
2518                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
2519                                             .AND.  ( nc_accr > eps_mr ) ) THEN
2520!
2521!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2522                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
2523!
2524!--                Universal function for accretion process (Seifert and
2525!--                Beheng, 2001):
2526                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2527
2528!
2529!--                Mean weight of cloud drops
2530                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
2531!
2532!--                Parameterized turbulence effects on autoconversion (Seifert,
2533!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2534!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2535                   IF ( collision_turbulence )  THEN
2536                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
2537                                       MIN( 600.0_wp,                          &
2538                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
2539                                     )
2540                   ELSE
2541                      k_cr = k_cr0
2542                   ENDIF
2543!
2544!--                Accretion rate (Seifert and Beheng, 2006):
2545                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
2546                          SQRT( rho_surface * hyrho(k) )
2547                   accr = MIN( accr, qc(k,j,i) / dt_micro )
2548
2549                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
2550                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
2551                   IF ( microphysics_morrison )  THEN
2552                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
2553                                         accr / xc * hyrho(k) * dt_micro * flag)
2554                   ENDIF
2555
2556                ENDIF
2557
2558             ENDDO
2559          ENDDO
2560       ENDDO
2561
2562       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
2563
2564    END SUBROUTINE accretion
2565
2566
2567!------------------------------------------------------------------------------!
2568! Description:
2569! ------------
2570!> Collisional breakup rate (Seifert, 2008).
2571!------------------------------------------------------------------------------!
2572    SUBROUTINE selfcollection_breakup
2573
2574       IMPLICIT NONE
2575
2576       INTEGER(iwp) ::  i                 !<
2577       INTEGER(iwp) ::  j                 !<
2578       INTEGER(iwp) ::  k                 !<
2579
2580       REAL(wp)     ::  breakup           !<
2581       REAL(wp)     ::  dr                !<
2582       REAL(wp)     ::  flag              !< flag to mask topography grid points
2583       REAL(wp)     ::  phi_br            !<
2584       REAL(wp)     ::  selfcoll          !<
2585
2586       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
2587
2588       DO  i = nxlg, nxrg
2589          DO  j = nysg, nyng
2590             DO  k = nzb+1, nzt
2591!
2592!--             Predetermine flag to mask topography
2593                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2594
2595                IF ( qr(k,j,i) > eps_sb )  THEN
2596!
2597!--                Selfcollection rate (Seifert and Beheng, 2001):
2598                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
2599                              SQRT( hyrho(k) * rho_surface )
2600!
2601!--                Weight averaged diameter of rain drops:
2602                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2603                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2604!
2605!--                Collisional breakup rate (Seifert, 2008):
2606                   IF ( dr >= 0.3E-3_wp )  THEN
2607                      phi_br  = k_br * ( dr - 1.1E-3_wp )
2608                      breakup = selfcoll * ( phi_br + 1.0_wp )
2609                   ELSE
2610                      breakup = 0.0_wp
2611                   ENDIF
2612
2613                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
2614                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
2615
2616                ENDIF
2617             ENDDO
2618          ENDDO
2619       ENDDO
2620
2621       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
2622
2623    END SUBROUTINE selfcollection_breakup
2624
2625
2626!------------------------------------------------------------------------------!
2627! Description:
2628! ------------
2629!> Evaporation of precipitable water. Condensation is neglected for
2630!> precipitable water.
2631!------------------------------------------------------------------------------!
2632    SUBROUTINE evaporation_rain
2633
2634       IMPLICIT NONE
2635
2636       INTEGER(iwp) ::  i                 !<
2637       INTEGER(iwp) ::  j                 !<
2638       INTEGER(iwp) ::  k                 !<
2639
2640       REAL(wp)     ::  dr                !<
2641       REAL(wp)     ::  evap              !<
2642       REAL(wp)     ::  evap_nr           !<
2643       REAL(wp)     ::  f_vent            !<
2644       REAL(wp)     ::  flag              !< flag to mask topography grid points
2645       REAL(wp)     ::  g_evap            !<
2646       REAL(wp)     ::  lambda_r          !<
2647       REAL(wp)     ::  mu_r              !<
2648       REAL(wp)     ::  mu_r_2            !<
2649       REAL(wp)     ::  mu_r_5d2          !<
2650       REAL(wp)     ::  nr_0              !<
2651       REAL(wp)     ::  temp              !<
2652       REAL(wp)     ::  xr                !<
2653
2654       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
2655
2656       DO  i = nxlg, nxrg
2657          DO  j = nysg, nyng
2658             DO  k = nzb+1, nzt
2659!
2660!--             Predetermine flag to mask topography
2661                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2662
2663                IF ( qr(k,j,i) > eps_sb )  THEN
2664
2665!
2666!--                Call calculation of supersaturation 
2667                   CALL supersaturation ( i, j, k )
2668!
2669!--                Evaporation needs only to be calculated in subsaturated regions
2670                   IF ( sat < 0.0_wp )  THEN
2671!
2672!--                   Actual temperature:
2673                      temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2674
2675                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
2676                                          l_v / ( thermal_conductivity_l * temp ) &
2677                                          + r_v * temp / ( diff_coeff_l * e_s )   &
2678                                        )
2679!
2680!--                   Mean weight of rain drops
2681                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
2682!
2683!--                   Weight averaged diameter of rain drops:
2684                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2685!
2686!--                   Compute ventilation factor and intercept parameter
2687!--                   (Seifert and Beheng, 2006; Seifert, 2008):
2688                      IF ( ventilation_effect )  THEN
2689!
2690!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
2691!--                      2005; Stevens and Seifert, 2008):
2692                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
2693                                                          ( dr - 1.4E-3_wp ) ) )
2694!
2695!--                      Slope parameter of gamma distribution (Seifert, 2008):
2696                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
2697                                      ( mu_r + 1.0_wp )                        &
2698                                    )**( 1.0_wp / 3.0_wp ) / dr
2699
2700                         mu_r_2   = mu_r + 2.0_wp
2701                         mu_r_5d2 = mu_r + 2.5_wp
2702
2703                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
2704                                  lambda_r**( -mu_r_2 ) + b_vent *              &
2705                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
2706                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
2707                                  ( 1.0_wp -                                    &
2708                                    0.5_wp * ( b_term / a_term ) *              &
2709                                    ( lambda_r / ( c_term + lambda_r )          &
2710                                    )**mu_r_5d2 -                               &
2711                                    0.125_wp * ( b_term / a_term )**2 *         &
2712                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
2713                                    )**mu_r_5d2 -                               &
2714                                    0.0625_wp * ( b_term / a_term )**3 *        &
2715                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
2716                                    )**mu_r_5d2 -                               &
2717                                    0.0390625_wp * ( b_term / a_term )**4 *     &
2718                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
2719                                    )**mu_r_5d2                                 &
2720                                  )
2721
2722                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
2723                                  gamm( mu_r + 1.0_wp )
2724                      ELSE
2725                         f_vent = 1.0_wp
2726                         nr_0   = nr(k,j,i) * dr
2727                      ENDIF
2728   !
2729   !--                Evaporation rate of rain water content (Seifert and
2730   !--                Beheng, 2006):
2731                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
2732                                hyrho(k)
2733                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
2734                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
2735                                     -nr(k,j,i) / dt_micro )
2736
2737                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
2738                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
2739
2740                   ENDIF
2741                ENDIF
2742
2743             ENDDO
2744          ENDDO
2745       ENDDO
2746
2747       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
2748
2749    END SUBROUTINE evaporation_rain
2750
2751
2752!------------------------------------------------------------------------------!
2753! Description:
2754! ------------
2755!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2756!------------------------------------------------------------------------------!
2757    SUBROUTINE sedimentation_cloud
2758
2759
2760       IMPLICIT NONE
2761
2762       INTEGER(iwp) ::  i             !<
2763       INTEGER(iwp) ::  j             !<
2764       INTEGER(iwp) ::  k             !<
2765
2766       REAL(wp) ::  flag              !< flag to mask topography grid points
2767       REAL(wp) ::  nc_sedi           !<
2768
2769       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
2770       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
2771
2772
2773       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
2774
2775       sed_qc(nzt+1) = 0.0_wp
2776       sed_nc(nzt+1) = 0.0_wp
2777
2778       DO  i = nxlg, nxrg
2779          DO  j = nysg, nyng
2780             DO  k = nzt, nzb+1, -1
2781!
2782!--             Predetermine flag to mask topography
2783                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2784
2785                IF ( microphysics_morrison ) THEN
2786                   nc_sedi = nc(k,j,i)
2787                ELSE
2788                   nc_sedi = nc_const
2789                ENDIF
2790
2791!
2792!--             Sedimentation fluxes for number concentration are only calculated
2793!--             for cloud_scheme = 'morrison'
2794                IF ( microphysics_morrison ) THEN
2795                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
2796                      sed_nc(k) = sed_qc_const *                               &
2797                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
2798                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
2799                   ELSE
2800                      sed_nc(k) = 0.0_wp
2801                   ENDIF
2802
2803                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
2804                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
2805                                  ) * flag
2806
2807                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
2808                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
2809                ENDIF
2810
2811                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
2812                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
2813                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
2814                                                                           flag
2815                ELSE
2816                   sed_qc(k) = 0.0_wp
2817                ENDIF
2818
2819                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
2820                                            dt_micro + sed_qc(k+1)             &
2821                               ) * flag
2822
2823                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
2824                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2825                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
2826                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2827                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
2828                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
2829                                        d_exner(k) * dt_micro            * flag
2830
2831!
2832!--             Compute the precipitation rate due to cloud (fog) droplets
2833                IF ( call_microphysics_at_all_substeps )  THEN
2834                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
2835                                * weight_substep(intermediate_timestep_count)  &
2836                                * flag
2837                ELSE
2838                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2839                ENDIF
2840
2841             ENDDO
2842          ENDDO
2843       ENDDO
2844
2845       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
2846
2847    END SUBROUTINE sedimentation_cloud
2848
2849
2850!------------------------------------------------------------------------------!
2851! Description:
2852! ------------
2853!> Computation of sedimentation flux. Implementation according to Stevens
2854!> and Seifert (2008). Code is based on UCLA-LES.
2855!------------------------------------------------------------------------------!
2856    SUBROUTINE sedimentation_rain
2857
2858       IMPLICIT NONE
2859
2860       INTEGER(iwp) ::  i             !< running index x direction
2861       INTEGER(iwp) ::  j             !< running index y direction
2862       INTEGER(iwp) ::  k             !< running index z direction
2863       INTEGER(iwp) ::  k_run         !<
2864       INTEGER(iwp) ::  m             !< running index surface elements
2865       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2866       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2867
2868       REAL(wp)     ::  c_run                      !<
2869       REAL(wp)     ::  d_max                      !<
2870       REAL(wp)     ::  d_mean                     !<
2871       REAL(wp)     ::  d_min                      !<
2872       REAL(wp)     ::  dr                         !<
2873       REAL(wp)     ::  flux                       !<
2874       REAL(wp)     ::  flag                       !< flag to mask topography grid points
2875       REAL(wp)     ::  lambda_r                   !<
2876       REAL(wp)     ::  mu_r                       !<
2877       REAL(wp)     ::  z_run                      !<
2878
2879       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2880       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2881       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2882       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2883       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2884       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2885       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2886       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
2887
2888       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
2889
2890!
2891!--    Compute velocities
2892       DO  i = nxlg, nxrg
2893          DO  j = nysg, nyng
2894             DO  k = nzb+1, nzt
2895!
2896!--             Predetermine flag to mask topography
2897                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2898
2899                IF ( qr(k,j,i) > eps_sb )  THEN
2900!
2901!--                Weight averaged diameter of rain drops:
2902                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2903                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2904!
2905!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2906!--                Stevens and Seifert, 2008):
2907                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
2908                                                     ( dr - 1.4E-3_wp ) ) )
2909!
2910!--                Slope parameter of gamma distribution (Seifert, 2008):
2911                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
2912                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2913
2914                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2915                                               a_term - b_term * ( 1.0_wp +    &
2916                                                  c_term /                     &
2917                                                  lambda_r )**( -1.0_wp *      &
2918                                                  ( mu_r + 1.0_wp ) )          &
2919                                              )                                &
2920                                ) * flag
2921
2922                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2923                                               a_term - b_term * ( 1.0_wp +    &
2924                                                  c_term /                     &
2925                                                  lambda_r )**( -1.0_wp *      &
2926                                                  ( mu_r + 4.0_wp ) )          &
2927                                             )                                 &
2928                                ) * flag
2929                ELSE
2930                   w_nr(k) = 0.0_wp
2931                   w_qr(k) = 0.0_wp
2932                ENDIF
2933             ENDDO
2934!
2935!--          Adjust boundary values using surface data type.
2936!--          Upward-facing
2937             surf_s = bc_h(0)%start_index(j,i)
2938             surf_e = bc_h(0)%end_index(j,i)
2939             DO  m = surf_s, surf_e
2940                k         = bc_h(0)%k(m)
2941                w_nr(k-1) = w_nr(k)
2942                w_qr(k-1) = w_qr(k)
2943             ENDDO
2944!
2945!--          Downward-facing
2946             surf_s = bc_h(1)%start_index(j,i)
2947             surf_e = bc_h(1)%end_index(j,i)
2948             DO  m = surf_s, surf_e
2949                k         = bc_h(1)%k(m)
2950                w_nr(k+1) = w_nr(k)
2951                w_qr(k+1) = w_qr(k)
2952             ENDDO
2953!
2954!--          Model top boundary value
2955             w_nr(nzt+1) = 0.0_wp
2956             w_qr(nzt+1) = 0.0_wp
2957!
2958!--          Compute Courant number
2959             DO  k = nzb+1, nzt
2960!
2961!--             Predetermine flag to mask topography
2962                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2963
2964                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
2965                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
2966                          dt_micro * ddzu(k) * flag
2967                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
2968                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
2969                          dt_micro * ddzu(k) * flag
2970             ENDDO
2971!
2972!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2973             IF ( limiter_sedimentation )  THEN
2974
2975                DO k = nzb+1, nzt
2976!
2977!--                Predetermine flag to mask topography
2978                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2979
2980                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
2981                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
2982                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
2983
2984                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
2985                                                              2.0_wp * d_max,  &
2986                                                              ABS( d_mean ) )  &
2987                                                      * flag
2988
2989                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
2990                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
2991                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
2992
2993                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
2994                                                              2.0_wp * d_max,  &
2995                                                              ABS( d_mean ) )
2996                ENDDO
2997
2998             ELSE
2999
3000                nr_slope = 0.0_wp
3001                qr_slope = 0.0_wp
3002
3003             ENDIF
3004
3005             sed_nr(nzt+1) = 0.0_wp
3006             sed_qr(nzt+1) = 0.0_wp
3007!
3008!--          Compute sedimentation flux
3009             DO  k = nzt, nzb+1, -1
3010!
3011!--             Predetermine flag to mask topography
3012                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3013!
3014!--             Sum up all rain drop number densities which contribute to the flux
3015!--             through k-1/2
3016                flux  = 0.0_wp
3017                z_run = 0.0_wp ! height above z(k)
3018                k_run = k
3019                c_run = MIN( 1.0_wp, c_nr(k) )
3020                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3021                   flux  = flux + hyrho(k_run) *                               &
3022                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
3023                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
3024                                              * flag
3025                   z_run = z_run + dzu(k_run) * flag
3026                   k_run = k_run + 1          * flag
3027                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
3028                                              * flag
3029                ENDDO
3030!
3031!--             It is not allowed to sediment more rain drop number density than
3032!--             available
3033                flux = MIN( flux,                                              &
3034                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
3035                            dt_micro                                           &
3036                          )
3037
3038                sed_nr(k) = flux / dt_micro * flag
3039                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
3040                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3041!
3042!--             Sum up all rain water content which contributes to the flux
3043!--             through k-1/2
3044                flux  = 0.0_wp
3045                z_run = 0.0_wp ! height above z(k)
3046                k_run = k
3047                c_run = MIN( 1.0_wp, c_qr(k) )
3048
3049                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3050
3051                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
3052                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
3053                                  0.5_wp ) * c_run * dzu(k_run) * flag
3054                   z_run = z_run + dzu(k_run)                   * flag
3055                   k_run = k_run + 1                            * flag
3056                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
3057                                                                * flag
3058
3059                ENDDO
3060!
3061!--             It is not allowed to sediment more rain water content than
3062!--             available
3063                flux = MIN( flux,                                              &
3064                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
3065                            dt_micro                                           &
3066                          )
3067
3068                sed_qr(k) = flux / dt_micro * flag
3069
3070                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
3071                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3072                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
3073                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3074                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
3075                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
3076                                        d_exner(k) * dt_micro            * flag
3077!
3078!--             Compute the rain rate
3079                IF ( call_microphysics_at_all_substeps )  THEN
3080                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
3081                                * weight_substep(intermediate_timestep_count) &
3082                                * flag
3083                ELSE
3084                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3085                ENDIF
3086
3087             ENDDO
3088          ENDDO
3089       ENDDO
3090
3091       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
3092
3093    END SUBROUTINE sedimentation_rain
3094
3095
3096!------------------------------------------------------------------------------!
3097! Description:
3098! ------------
3099!> Computation of the precipitation amount due to gravitational settling of
3100!> rain and cloud (fog) droplets
3101!------------------------------------------------------------------------------!
3102    SUBROUTINE calc_precipitation_amount
3103
3104       IMPLICIT NONE
3105
3106       INTEGER(iwp) ::  i             !< running index x direction
3107       INTEGER(iwp) ::  j             !< running index y direction
3108       INTEGER(iwp) ::  k             !< running index y direction
3109       INTEGER(iwp) ::  m             !< running index surface elements
3110
3111       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3112            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
3113            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3114       THEN
3115!
3116!--       Run over all upward-facing surface elements, i.e. non-natural,
3117!--       natural and urban
3118          DO  m = 1, bc_h(0)%ns
3119             i = bc_h(0)%i(m)
3120             j = bc_h(0)%j(m)
3121             k = bc_h(0)%k(m)
3122             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
3123                                               prr(k,j,i) * hyrho(k) * dt_3d
3124          ENDDO
3125
3126       ENDIF
3127
3128    END SUBROUTINE calc_precipitation_amount
3129
3130
3131!------------------------------------------------------------------------------!
3132! Description:
3133! ------------
3134!> Control of microphysics for grid points i,j
3135!------------------------------------------------------------------------------!
3136
3137    SUBROUTINE bcm_actions_ij( i, j )
3138
3139       IMPLICIT NONE
3140
3141       INTEGER(iwp) ::  i                 !<
3142       INTEGER(iwp) ::  j                 !<
3143
3144       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
3145!
3146!--       Calculate vertical profile of the hydrostatic pressure (hyp)
3147          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
3148          d_exner = exner_function_invers(hyp)
3149          exner = 1.0_wp / exner_function_invers(hyp)
3150          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
3151!
3152!--       Compute reference density
3153          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
3154       ENDIF
3155
3156!
3157!--    Compute length of time step
3158       IF ( call_microphysics_at_all_substeps )  THEN
3159          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
3160       ELSE
3161          dt_micro = dt_3d
3162       ENDIF
3163!
3164!--    Reset precipitation rate
3165       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
3166
3167!
3168!--    Compute cloud physics
3169       IF( microphysics_kessler )  THEN
3170
3171          CALL autoconversion_kessler( i,j )
3172          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3173
3174       ELSEIF ( microphysics_seifert )  THEN
3175
3176          CALL adjust_cloud( i,j )
3177          IF ( microphysics_morrison ) CALL activation( i,j )
3178          IF ( microphysics_morrison ) CALL condensation( i,j )
3179          CALL autoconversion( i,j )
3180          CALL accretion( i,j )
3181          CALL selfcollection_breakup( i,j )
3182          CALL evaporation_rain( i,j )
3183          CALL sedimentation_rain( i,j )
3184          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3185
3186       ENDIF
3187
3188       CALL calc_precipitation_amount( i,j )
3189
3190    END SUBROUTINE bcm_actions_ij
3191
3192!------------------------------------------------------------------------------!
3193! Description:
3194! ------------
3195!> Adjust number of raindrops to avoid nonlinear effects in
3196!> sedimentation and evaporation of rain drops due to too small or
3197!> too big weights of rain drops (Stevens and Seifert, 2008).
3198!> The same procedure is applied to cloud droplets if they are determined
3199!> prognostically. Call for grid point i,j
3200!------------------------------------------------------------------------------!
3201    SUBROUTINE adjust_cloud_ij( i, j )
3202
3203       IMPLICIT NONE
3204
3205       INTEGER(iwp) ::  i                 !<
3206       INTEGER(iwp) ::  j                 !<
3207       INTEGER(iwp) ::  k                 !<
3208
3209       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
3210
3211       DO  k = nzb+1, nzt
3212!
3213!--       Predetermine flag to mask topography
3214          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3215
3216          IF ( qr(k,j,i) <= eps_sb )  THEN
3217             qr(k,j,i) = 0.0_wp
3218             nr(k,j,i) = 0.0_wp
3219          ELSE
3220!
3221!--          Adjust number of raindrops to avoid nonlinear effects in
3222!--          sedimentation and evaporation of rain drops due to too small or
3223!--          too big weights of rain drops (Stevens and Seifert, 2008).
3224             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
3225                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
3226             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
3227                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
3228             ENDIF
3229
3230          ENDIF
3231
3232          IF ( microphysics_morrison ) THEN
3233             IF ( qc(k,j,i) <= eps_sb )  THEN
3234                qc(k,j,i) = 0.0_wp
3235                nc(k,j,i) = 0.0_wp
3236             ELSE
3237                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
3238                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
3239                ENDIF
3240             ENDIF
3241          ENDIF
3242
3243       ENDDO
3244
3245    END SUBROUTINE adjust_cloud_ij
3246
3247!------------------------------------------------------------------------------!
3248! Description:
3249! ------------
3250!> Calculate number of activated condensation nucleii after simple activation
3251!> scheme of Twomey, 1959.
3252!------------------------------------------------------------------------------!
3253    SUBROUTINE activation_ij( i, j )
3254
3255       IMPLICIT NONE
3256
3257       INTEGER(iwp) ::  i                 !<
3258       INTEGER(iwp) ::  j                 !<
3259       INTEGER(iwp) ::  k                 !<
3260
3261       REAL(wp)     ::  activ             !<
3262       REAL(wp)     ::  afactor           !<
3263       REAL(wp)     ::  beta_act          !<
3264       REAL(wp)     ::  bfactor           !<
3265       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3266       REAL(wp)     ::  k_act             !<
3267       REAL(wp)     ::  n_act             !<
3268       REAL(wp)     ::  n_ccn             !<
3269       REAL(wp)     ::  s_0               !<
3270       REAL(wp)     ::  sat_max           !<
3271       REAL(wp)     ::  sigma             !<
3272       REAL(wp)     ::  sigma_act         !<
3273
3274       DO  k = nzb+1, nzt
3275!
3276!--       Predetermine flag to mask topography
3277          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3278!
3279!--       Call calculation of supersaturation 
3280          CALL supersaturation ( i, j, k )
3281!
3282!--       Prescribe parameters for activation
3283!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
3284          k_act  = 0.7_wp
3285          activ  = 0.0_wp
3286
3287          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
3288!
3289!--          Compute the number of activated Aerosols
3290!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
3291             n_act     = na_init * sat**k_act
3292!
3293!--          Compute the number of cloud droplets
3294!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3295!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
3296
3297!
3298!--          Compute activation rate after Khairoutdinov and Kogan
3299!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
3300             sat_max = 0.8_wp / 100.0_wp
3301             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
3302                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
3303                  dt_micro
3304
3305             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
3306          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
3307!
3308!--          Curvature effect (afactor) with surface tension
3309!--          parameterization by Straka (2009)
3310             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
3311             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
3312!
3313!--          Solute effect (bfactor)
3314             bfactor = vanthoff * molecular_weight_of_water *                  &
3315                  rho_s / ( molecular_weight_of_solute * rho_l )
3316
3317!
3318!--          Prescribe power index that describes the soluble fraction
3319!--          of an aerosol particle (beta).
3320!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3321             beta_act  = 0.5_wp
3322             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
3323!
3324!--          Calculate mean geometric supersaturation (s_0) with
3325!--          parameterization by Khvorostyanov and Curry (2006)
3326             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
3327               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
3328
3329!
3330!--          Calculate number of activated CCN as a function of
3331!--          supersaturation and taking Koehler theory into account
3332!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
3333             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
3334                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
3335             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
3336
3337             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
3338          ENDIF
3339
3340       ENDDO
3341
3342    END SUBROUTINE activation_ij
3343
3344!------------------------------------------------------------------------------!
3345! Description:
3346! ------------
3347!> Calculate condensation rate for cloud water content (after Khairoutdinov and
3348!> Kogan, 2000).
3349!------------------------------------------------------------------------------!
3350    SUBROUTINE condensation_ij( i, j )
3351
3352       IMPLICIT NONE
3353
3354       INTEGER(iwp) ::  i                 !<
3355       INTEGER(iwp) ::  j                 !<
3356       INTEGER(iwp) ::  k                 !<
3357
3358       REAL(wp)     ::  cond              !<
3359       REAL(wp)     ::  cond_max          !<
3360       REAL(wp)     ::  dc                !<
3361       REAL(wp)     ::  evap              !<
3362       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3363       REAL(wp)     ::  g_fac             !<
3364       REAL(wp)     ::  nc_0              !<
3365       REAL(wp)     ::  temp              !<
3366       REAL(wp)     ::  xc                !<
3367
3368
3369       DO  k = nzb+1, nzt
3370!
3371!--       Predetermine flag to mask topography
3372          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3373!
3374!--       Call calculation of supersaturation 
3375          CALL supersaturation ( i, j, k )
3376!
3377!--       Actual temperature:
3378          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3379
3380          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
3381                              l_v / ( thermal_conductivity_l * temp )    &
3382                              + r_v * temp / ( diff_coeff_l * e_s )      &
3383                            )
3384!
3385!--       Mean weight of cloud drops
3386          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
3387          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
3388!
3389!--       Weight averaged diameter of cloud drops:
3390          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3391!
3392!--       Integral diameter of cloud drops
3393          nc_0 = nc(k,j,i) * dc
3394!
3395!--       Condensation needs only to be calculated in supersaturated regions
3396          IF ( sat > 0.0_wp )  THEN
3397!
3398!--          Condensation rate of cloud water content
3399!--          after KK scheme.
3400!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
3401             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3402             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
3403             cond      = MIN( cond, cond_max / dt_micro )
3404
3405             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
3406          ELSEIF ( sat < 0.0_wp ) THEN
3407             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3408             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
3409
3410             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
3411          ENDIF
3412       ENDDO
3413
3414    END SUBROUTINE condensation_ij
3415
3416
3417!------------------------------------------------------------------------------!
3418! Description:
3419! ------------
3420!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
3421!------------------------------------------------------------------------------!
3422    SUBROUTINE autoconversion_ij( i, j )
3423
3424       IMPLICIT NONE
3425
3426       INTEGER(iwp) ::  i                 !<
3427       INTEGER(iwp) ::  j                 !<
3428       INTEGER(iwp) ::  k                 !<
3429
3430       REAL(wp)     ::  alpha_cc          !<
3431       REAL(wp)     ::  autocon           !<
3432       REAL(wp)     ::  dissipation       !<
3433       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3434       REAL(wp)     ::  k_au              !<
3435       REAL(wp)     ::  l_mix             !<
3436       REAL(wp)     ::  nc_auto           !<
3437       REAL(wp)     ::  nu_c              !<
3438       REAL(wp)     ::  phi_au            !<
3439       REAL(wp)     ::  r_cc              !<
3440       REAL(wp)     ::  rc                !<
3441       REAL(wp)     ::  re_lambda         !<
3442       REAL(wp)     ::  sigma_cc          !<
3443       REAL(wp)     ::  tau_cloud         !<
3444       REAL(wp)     ::  xc                !<
3445
3446       DO  k = nzb+1, nzt
3447!
3448!--       Predetermine flag to mask topography
3449          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3450          IF ( microphysics_morrison ) THEN
3451             nc_auto = nc(k,j,i)
3452          ELSE
3453             nc_auto = nc_const
3454          ENDIF
3455
3456          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
3457
3458             k_au = k_cc / ( 20.0_wp * x0 )
3459!
3460!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3461!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
3462             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
3463                              0.0_wp )
3464!
3465!--          Universal function for autoconversion process
3466!--          (Seifert and Beheng, 2006):
3467             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
3468!
3469!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
3470!--          (Use constant nu_c = 1.0_wp instead?)
3471             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
3472!
3473!--          Mean weight of cloud droplets:
3474             xc = hyrho(k) * qc(k,j,i) / nc_auto
3475!
3476!--          Parameterized turbulence effects on autoconversion (Seifert,
3477!--          Nuijens and Stevens, 2010)
3478             IF ( collision_turbulence )  THEN
3479!
3480!--             Weight averaged radius of cloud droplets:
3481                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3482
3483                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
3484                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
3485                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
3486!
3487!--             Mixing length (neglecting distance to ground and stratification)
3488                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
3489!
3490!--             Limit dissipation rate according to Seifert, Nuijens and
3491!--             Stevens (2010)
3492                dissipation = MIN( 0.06_wp, diss(k,j,i) )
3493!
3494!--             Compute Taylor-microscale Reynolds number:
3495                re_lambda = 6.0_wp / 11.0_wp *                                 &
3496                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
3497                            SQRT( 15.0_wp / kin_vis_air ) *                    &
3498                            dissipation**( 1.0_wp / 6.0_wp )
3499!
3500!--             The factor of 1.0E4 is needed to convert the dissipation rate
3501!--             from m2 s-3 to cm2 s-3.
3502                k_au = k_au * ( 1.0_wp +                                       &
3503                                dissipation * 1.0E4_wp *                       &
3504                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
3505                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
3506                                                  sigma_cc )**2                &
3507                                                ) + beta_cc                    &
3508                                )                                              &
3509                              )
3510             ENDIF
3511!
3512!--          Autoconversion rate (Seifert and Beheng, 2006):
3513             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
3514                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
3515                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
3516                       rho_surface
3517             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
3518
3519             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
3520             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
3521             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
3522             IF ( microphysics_morrison ) THEN
3523                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
3524                                    autocon / x0 * hyrho(k) * dt_micro * flag )
3525             ENDIF
3526
3527          ENDIF
3528
3529       ENDDO
3530
3531    END SUBROUTINE autoconversion_ij
3532
3533!------------------------------------------------------------------------------!
3534! Description:
3535! ------------
3536!> Autoconversion process (Kessler, 1969).
3537!------------------------------------------------------------------------------!
3538    SUBROUTINE autoconversion_kessler_ij( i, j )
3539
3540
3541       IMPLICIT NONE
3542
3543       INTEGER(iwp) ::  i      !<
3544       INTEGER(iwp) ::  j      !<
3545       INTEGER(iwp) ::  k      !<
3546       INTEGER(iwp) ::  k_wall !< topography top index
3547
3548       REAL(wp)    ::  dqdt_precip !<
3549       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
3550
3551!
3552!--    Determine vertical index of topography top
3553       k_wall = get_topography_top_index_ji( j, i, 's' )
3554       DO  k = nzb+1, nzt
3555!
3556!--       Predetermine flag to mask topography
3557          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3558
3559          IF ( qc(k,j,i) > ql_crit )  THEN
3560             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
3561          ELSE
3562             dqdt_precip = 0.0_wp
3563          ENDIF
3564
3565          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
3566          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
3567          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
3568                                                                 * flag
3569
3570!
3571!--       Compute the rain rate (stored on surface grid point)
3572          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
3573
3574       ENDDO
3575
3576    END SUBROUTINE autoconversion_kessler_ij
3577
3578!------------------------------------------------------------------------------!
3579! Description:
3580! ------------
3581!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
3582!------------------------------------------------------------------------------!
3583    SUBROUTINE accretion_ij( i, j )
3584
3585       IMPLICIT NONE
3586
3587       INTEGER(iwp) ::  i                 !<
3588       INTEGER(iwp) ::  j                 !<
3589       INTEGER(iwp) ::  k                 !<
3590
3591       REAL(wp)     ::  accr              !<
3592       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3593       REAL(wp)     ::  k_cr              !<
3594       REAL(wp)     ::  nc_accr           !<
3595       REAL(wp)     ::  phi_ac            !<
3596       REAL(wp)     ::  tau_cloud         !<
3597       REAL(wp)     ::  xc                !<
3598
3599
3600       DO  k = nzb+1, nzt
3601!
3602!--       Predetermine flag to mask topography
3603          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3604          IF ( microphysics_morrison ) THEN
3605             nc_accr = nc(k,j,i)
3606          ELSE
3607             nc_accr = nc_const
3608          ENDIF
3609
3610          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
3611               ( nc_accr > eps_mr ) )  THEN
3612!
3613!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3614             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
3615!
3616!--          Universal function for accretion process
3617!--          (Seifert and Beheng, 2001):
3618             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
3619
3620!
3621!--          Mean weight of cloud drops
3622             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
3623!
3624!--          Parameterized turbulence effects on autoconversion (Seifert,
3625!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
3626!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
3627             IF ( collision_turbulence )  THEN
3628                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
3629                                 MIN( 600.0_wp,                                &
3630                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
3631                               )
3632             ELSE
3633                k_cr = k_cr0
3634             ENDIF
3635!
3636!--          Accretion rate (Seifert and Beheng, 2006):
3637             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
3638                    SQRT( rho_surface * hyrho(k) )
3639             accr = MIN( accr, qc(k,j,i) / dt_micro )
3640
3641             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
3642             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
3643             IF ( microphysics_morrison )  THEN
3644                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
3645                                             hyrho(k) * dt_micro * flag        &
3646                                         )
3647             ENDIF
3648
3649
3650          ENDIF
3651
3652       ENDDO
3653
3654    END SUBROUTINE accretion_ij
3655
3656
3657!------------------------------------------------------------------------------!
3658! Description:
3659! ------------
3660!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
3661!------------------------------------------------------------------------------!
3662    SUBROUTINE selfcollection_breakup_ij( i, j )
3663
3664       IMPLICIT NONE
3665
3666       INTEGER(iwp) ::  i                 !<
3667       INTEGER(iwp) ::  j                 !<
3668       INTEGER(iwp) ::  k                 !<
3669
3670       REAL(wp)     ::  breakup           !<
3671       REAL(wp)     ::  dr                !<
3672       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3673       REAL(wp)     ::  phi_br            !<
3674       REAL(wp)     ::  selfcoll          !<
3675
3676       DO  k = nzb+1, nzt
3677!
3678!--       Predetermine flag to mask topography
3679          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3680
3681          IF ( qr(k,j,i) > eps_sb )  THEN
3682!
3683!--          Selfcollection rate (Seifert and Beheng, 2001):
3684             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
3685!
3686!--          Weight averaged diameter of rain drops:
3687             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3688!
3689!--          Collisional breakup rate (Seifert, 2008):
3690             IF ( dr >= 0.3E-3_wp )  THEN
3691                phi_br  = k_br * ( dr - 1.1E-3_wp )
3692                breakup = selfcoll * ( phi_br + 1.0_wp )
3693             ELSE
3694                breakup = 0.0_wp
3695             ENDIF
3696
3697             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
3698             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
3699
3700          ENDIF
3701       ENDDO
3702
3703    END SUBROUTINE selfcollection_breakup_ij
3704
3705
3706!------------------------------------------------------------------------------!
3707! Description:
3708! ------------
3709!> Evaporation of precipitable water. Condensation is neglected for
3710!> precipitable water. Call for grid point i,j
3711!------------------------------------------------------------------------------!
3712    SUBROUTINE evaporation_rain_ij( i, j )
3713
3714       IMPLICIT NONE
3715
3716       INTEGER(iwp) ::  i                 !<
3717       INTEGER(iwp) ::  j                 !<
3718       INTEGER(iwp) ::  k                 !<
3719
3720       REAL(wp)     ::  dr                !<
3721       REAL(wp)     ::  evap              !<
3722       REAL(wp)     ::  evap_nr           !<
3723       REAL(wp)     ::  f_vent            !<
3724       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3725       REAL(wp)     ::  g_evap            !<
3726       REAL(wp)     ::  lambda_r          !<
3727       REAL(wp)     ::  mu_r              !<
3728       REAL(wp)     ::  mu_r_2            !<
3729       REAL(wp)     ::  mu_r_5d2          !<
3730       REAL(wp)     ::  nr_0              !<
3731       REAL(wp)     ::  temp              !<
3732       REAL(wp)     ::  xr                !<
3733
3734       DO  k = nzb+1, nzt
3735!
3736!--       Predetermine flag to mask topography
3737          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3738
3739          IF ( qr(k,j,i) > eps_sb )  THEN
3740!
3741!--          Call calculation of supersaturation 
3742             CALL supersaturation ( i, j, k )
3743!
3744!--          Evaporation needs only to be calculated in subsaturated regions
3745             IF ( sat < 0.0_wp )  THEN
3746!
3747!--             Actual temperature:
3748                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3749
3750                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
3751                                    ( thermal_conductivity_l * temp ) +        &
3752                                    r_v * temp / ( diff_coeff_l * e_s )        &
3753                                  )
3754!
3755!--             Mean weight of rain drops
3756                xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
3757!
3758!--             Weight averaged diameter of rain drops:
3759                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
3760!
3761!--             Compute ventilation factor and intercept parameter
3762!--             (Seifert and Beheng, 2006; Seifert, 2008):
3763                IF ( ventilation_effect )  THEN
3764!
3765!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3766!--                Stevens and Seifert, 2008):
3767                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3768!
3769!--                Slope parameter of gamma distribution (Seifert, 2008):
3770                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
3771                                ( mu_r + 1.0_wp )                              &
3772                              )**( 1.0_wp / 3.0_wp ) / dr
3773
3774                   mu_r_2   = mu_r + 2.0_wp
3775                   mu_r_5d2 = mu_r + 2.5_wp
3776
3777                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
3778                            b_vent * schmidt_p_1d3 *                           &
3779                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
3780                            lambda_r**( -mu_r_5d2 ) *                          &
3781                            ( 1.0_wp -                                         &
3782                              0.5_wp * ( b_term / a_term ) *                   &
3783                              ( lambda_r / ( c_term + lambda_r )               &
3784                              )**mu_r_5d2 -                                    &
3785                              0.125_wp * ( b_term / a_term )**2 *              &
3786                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
3787                              )**mu_r_5d2 -                                    &
3788                              0.0625_wp * ( b_term / a_term )**3 *             &
3789                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
3790                              )**mu_r_5d2 -                                    &
3791                              0.0390625_wp * ( b_term / a_term )**4 *          &
3792                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
3793                              )**mu_r_5d2                                      &
3794                            )
3795
3796                   nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /           &
3797                            gamm( mu_r + 1.0_wp )
3798                ELSE
3799                   f_vent = 1.0_wp
3800                   nr_0   = nr(k,j,i) * dr
3801                ENDIF
3802!
3803!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
3804                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
3805                evap    = MAX( evap, -qr(k,j,i) / dt_micro )
3806                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
3807                               -nr(k,j,i) / dt_micro )
3808
3809                qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
3810                nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
3811
3812             ENDIF
3813          ENDIF
3814
3815       ENDDO
3816
3817    END SUBROUTINE evaporation_rain_ij
3818
3819
3820!------------------------------------------------------------------------------!
3821! Description:
3822! ------------
3823!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
3824!> Call for grid point i,j
3825!------------------------------------------------------------------------------!
3826    SUBROUTINE sedimentation_cloud_ij( i, j )
3827
3828       IMPLICIT NONE
3829
3830       INTEGER(iwp) ::  i             !<
3831       INTEGER(iwp) ::  j             !<
3832       INTEGER(iwp) ::  k             !<
3833
3834       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
3835       REAL(wp)     ::  nc_sedi !<
3836
3837       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
3838       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
3839
3840       sed_qc(nzt+1) = 0.0_wp
3841       sed_nc(nzt+1) = 0.0_wp
3842
3843
3844       DO  k = nzt, nzb+1, -1
3845!
3846!--       Predetermine flag to mask topography
3847          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3848          IF ( microphysics_morrison ) THEN
3849             nc_sedi = nc(k,j,i)
3850          ELSE
3851             nc_sedi = nc_const
3852          ENDIF
3853!
3854!--       Sedimentation fluxes for number concentration are only calculated
3855!--       for cloud_scheme = 'morrison'
3856          IF ( microphysics_morrison ) THEN
3857             IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
3858                sed_nc(k) = sed_qc_const *                                     &
3859                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
3860                            ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
3861             ELSE
3862                sed_nc(k) = 0.0_wp
3863             ENDIF
3864
3865             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
3866                              nc(k,j,i) / dt_micro + sed_nc(k+1)                &
3867                            ) * flag
3868
3869             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *               &
3870                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
3871          ENDIF
3872
3873          IF ( qc(k,j,i) > eps_sb  .AND.  nc_sedi > eps_mr )  THEN
3874             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
3875                         ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
3876          ELSE
3877             sed_qc(k) = 0.0_wp
3878          ENDIF
3879
3880          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /          &
3881                                      dt_micro + sed_qc(k+1)                   &
3882                         ) * flag
3883
3884          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3885                                hyrho(k) * dt_micro * flag
3886          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3887                                hyrho(k) * dt_micro * flag
3888          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3889                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
3890                                                                * flag
3891
3892!
3893!--       Compute the precipitation rate of cloud (fog) droplets
3894          IF ( call_microphysics_at_all_substeps )  THEN
3895             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
3896                              weight_substep(intermediate_timestep_count) * flag
3897          ELSE
3898             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
3899          ENDIF
3900
3901       ENDDO
3902
3903    END SUBROUTINE sedimentation_cloud_ij
3904
3905
3906!------------------------------------------------------------------------------!
3907! Description:
3908! ------------
3909!> Computation of sedimentation flux. Implementation according to Stevens
3910!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
3911!------------------------------------------------------------------------------!
3912    SUBROUTINE sedimentation_rain_ij( i, j )
3913
3914       IMPLICIT NONE
3915
3916       INTEGER(iwp) ::  i             !< running index x direction
3917       INTEGER(iwp) ::  j             !< running index y direction
3918       INTEGER(iwp) ::  k             !< running index z direction
3919       INTEGER(iwp) ::  k_run         !<
3920       INTEGER(iwp) ::  m             !< running index surface elements
3921       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
3922       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
3923
3924       REAL(wp)     ::  c_run                      !<
3925       REAL(wp)     ::  d_max                      !<
3926       REAL(wp)     ::  d_mean                     !<
3927       REAL(wp)     ::  d_min                      !<
3928       REAL(wp)     ::  dr                         !<
3929       REAL(wp)     ::  flux                       !<
3930       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
3931       REAL(wp)     ::  lambda_r                   !<
3932       REAL(wp)     ::  mu_r                       !<
3933       REAL(wp)     ::  z_run                      !<
3934
3935       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
3936       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
3937       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
3938       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
3939       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
3940       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
3941       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
3942       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
3943
3944!
3945!--    Compute velocities
3946       DO  k = nzb+1, nzt
3947!
3948!--       Predetermine flag to mask topography
3949          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3950
3951          IF ( qr(k,j,i) > eps_sb )  THEN
3952!
3953!--          Weight averaged diameter of rain drops:
3954             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3955!
3956!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3957!--          Stevens and Seifert, 2008):
3958             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3959!
3960!--          Slope parameter of gamma distribution (Seifert, 2008):
3961             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
3962                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
3963
3964             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
3965                                         a_term - b_term * ( 1.0_wp +          &
3966                                            c_term / lambda_r )**( -1.0_wp *   &
3967                                            ( mu_r + 1.0_wp ) )                &
3968                                        )                                      &
3969                          ) * flag
3970             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
3971                                         a_term - b_term * ( 1.0_wp +          &
3972                                            c_term / lambda_r )**( -1.0_wp *   &
3973                                            ( mu_r + 4.0_wp ) )                &
3974                                       )                                       &
3975                          ) * flag
3976          ELSE
3977             w_nr(k) = 0.0_wp
3978             w_qr(k) = 0.0_wp
3979          ENDIF
3980       ENDDO
3981!
3982!--    Adjust boundary values using surface data type.
3983!--    Upward facing non-natural
3984       surf_s = bc_h(0)%start_index(j,i)
3985       surf_e = bc_h(0)%end_index(j,i)
3986       DO  m = surf_s, surf_e
3987          k         = bc_h(0)%k(m)
3988          w_nr(k-1) = w_nr(k)
3989          w_qr(k-1) = w_qr(k)
3990       ENDDO
3991!
3992!--    Downward facing non-natural
3993       surf_s = bc_h(1)%start_index(j,i)
3994       surf_e = bc_h(1)%end_index(j,i)
3995       DO  m = surf_s, surf_e
3996          k         = bc_h(1)%k(m)
3997          w_nr(k+1) = w_nr(k)
3998          w_qr(k+1) = w_qr(k)
3999       ENDDO
4000!
4001!--    Neumann boundary condition at model top
4002       w_nr(nzt+1) = 0.0_wp
4003       w_qr(nzt+1) = 0.0_wp
4004!
4005!--    Compute Courant number
4006       DO  k = nzb+1, nzt
4007!
4008!--       Predetermine flag to mask topography
4009          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4010
4011          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
4012                    dt_micro * ddzu(k) * flag
4013          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
4014                    dt_micro * ddzu(k) * flag
4015       ENDDO
4016!
4017!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
4018       IF ( limiter_sedimentation )  THEN
4019
4020          DO k = nzb+1, nzt
4021!
4022!--          Predetermine flag to mask topography
4023             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4024
4025             d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
4026             d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
4027             d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
4028
4029             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4030                                                        2.0_wp * d_max,        &
4031                                                        ABS( d_mean ) ) * flag
4032
4033             d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
4034             d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
4035             d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
4036
4037             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4038                                                        2.0_wp * d_max,        &
4039                                                        ABS( d_mean ) ) * flag
4040          ENDDO
4041
4042       ELSE
4043
4044          nr_slope = 0.0_wp
4045          qr_slope = 0.0_wp
4046
4047       ENDIF
4048
4049       sed_nr(nzt+1) = 0.0_wp
4050       sed_qr(nzt+1) = 0.0_wp
4051!
4052!--    Compute sedimentation flux
4053       DO  k = nzt, nzb+1, -1
4054!
4055!--       Predetermine flag to mask topography
4056          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4057!
4058!--       Sum up all rain drop number densities which contribute to the flux
4059!--       through k-1/2
4060          flux  = 0.0_wp
4061          z_run = 0.0_wp ! height above z(k)
4062          k_run = k
4063          c_run = MIN( 1.0_wp, c_nr(k) )
4064          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4065             flux  = flux + hyrho(k_run) *                                     &
4066                     ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4067                     0.5_wp ) * c_run * dzu(k_run) * flag
4068             z_run = z_run + dzu(k_run)            * flag
4069             k_run = k_run + 1                     * flag
4070             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
4071          ENDDO
4072!
4073!--       It is not allowed to sediment more rain drop number density than
4074!--       available
4075          flux = MIN( flux,                                                    &
4076                      hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * dt_micro )
4077
4078          sed_nr(k) = flux / dt_micro * flag
4079          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
4080                                    hyrho(k) * dt_micro * flag
4081!
4082!--       Sum up all rain water content which contributes to the flux
4083!--       through k-1/2
4084          flux  = 0.0_wp
4085          z_run = 0.0_wp ! height above z(k)
4086          k_run = k
4087          c_run = MIN( 1.0_wp, c_qr(k) )
4088
4089          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4090
4091             flux  = flux + hyrho(k_run) *                                     &
4092                     ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4093                     0.5_wp ) * c_run * dzu(k_run) * flag
4094             z_run = z_run + dzu(k_run)            * flag
4095             k_run = k_run + 1                     * flag
4096             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
4097
4098          ENDDO
4099!
4100!--       It is not allowed to sediment more rain water content than available
4101          flux = MIN( flux,                                                    &
4102                      hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_micro )
4103
4104          sed_qr(k) = flux / dt_micro * flag
4105
4106          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4107                                hyrho(k) * dt_micro * flag
4108          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4109                                hyrho(k) * dt_micro * flag
4110          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4111                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
4112                                                                * flag
4113!
4114!--       Compute the rain rate
4115          IF ( call_microphysics_at_all_substeps )  THEN
4116             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
4117                          * weight_substep(intermediate_timestep_count) * flag
4118          ELSE
4119             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
4120          ENDIF
4121
4122       ENDDO
4123
4124    END SUBROUTINE sedimentation_rain_ij
4125
4126
4127!------------------------------------------------------------------------------!
4128! Description:
4129! ------------
4130!> This subroutine computes the precipitation amount due to gravitational
4131!> settling of rain and cloud (fog) droplets
4132!------------------------------------------------------------------------------!
4133    SUBROUTINE calc_precipitation_amount_ij( i, j )
4134
4135       IMPLICIT NONE
4136
4137       INTEGER(iwp) ::  i             !< running index x direction
4138       INTEGER(iwp) ::  j             !< running index y direction
4139       INTEGER(iwp) ::  k             !< running index z direction
4140       INTEGER(iwp) ::  m             !< running index surface elements
4141       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
4142       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
4143
4144       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
4145            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
4146            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
4147       THEN
4148
4149          surf_s = bc_h(0)%start_index(j,i)
4150          surf_e = bc_h(0)%end_index(j,i)
4151          DO  m = surf_s, surf_e
4152             k                         = bc_h(0)%k(m)
4153             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
4154                                               prr(k,j,i) * hyrho(k) * dt_3d
4155          ENDDO
4156
4157       ENDIF
4158
4159    END SUBROUTINE calc_precipitation_amount_ij
4160
4161
4162!------------------------------------------------------------------------------!
4163! Description:
4164! ------------
4165!> Computation of the diagnostic supersaturation sat, actual liquid water
4166!< temperature t_l and saturation water vapor mixing ratio q_s
4167!------------------------------------------------------------------------------!
4168    SUBROUTINE supersaturation ( i,j,k )
4169
4170       IMPLICIT NONE
4171
4172       INTEGER(iwp) ::  i   !< running index
4173       INTEGER(iwp) ::  j   !< running index
4174       INTEGER(iwp) ::  k   !< running index
4175
4176       REAL(wp) ::  alpha   !< correction factor
4177!
4178!--    Actual liquid water temperature:
4179       t_l = exner(k) * pt(k,j,i)
4180!
4181!--    Calculate water vapor saturation pressure
4182       e_s = magnus( t_l )
4183!
4184!--    Computation of saturation mixing ratio:
4185       q_s   = rd_d_rv * e_s / ( hyp(k) - e_s )
4186!
4187!--    Correction factor
4188       alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l )
4189!
4190!--    Correction of the approximated value
4191!--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
4192       q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
4193!
4194!--    Supersaturation:
4195!--    Not in case of microphysics_kessler or microphysics_sat_adjust
4196!--    since qr is unallocated
4197       IF ( .NOT. microphysics_kessler  .AND.                                  &
4198            .NOT.  microphysics_sat_adjust )  THEN
4199          sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
4200       ENDIF
4201
4202    END SUBROUTINE supersaturation
4203
4204
4205!------------------------------------------------------------------------------!
4206! Description:
4207! ------------
4208!> Calculation of the liquid water content (0%-or-100%-scheme). This scheme is
4209!> used by the one and the two moment cloud physics scheme. Using the two moment
4210!> scheme, this calculation results in the cloud water content.
4211!------------------------------------------------------------------------------!
4212    SUBROUTINE calc_liquid_water_content
4213
4214
4215
4216       IMPLICIT NONE
4217
4218       INTEGER(iwp) ::  i !<
4219       INTEGER(iwp) ::  j !<
4220       INTEGER(iwp) ::  k !<
4221
4222       REAL(wp)     ::  flag !< flag to indicate first grid level above surface
4223
4224
4225       DO  i = nxlg, nxrg
4226          DO  j = nysg, nyng
4227             DO  k = nzb+1, nzt
4228!
4229!--             Predetermine flag to mask topography
4230                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4231
4232   !
4233   !--          Call calculation of supersaturation located
4234                CALL supersaturation( i, j, k )
4235
4236   !
4237   !--          Compute the liquid water content
4238                IF ( microphysics_seifert  .AND. .NOT. microphysics_morrison ) &
4239                THEN
4240                   IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
4241                      qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
4242                      ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
4243                   ELSE
4244                      IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
4245                      qc(k,j,i) = 0.0_wp
4246                      ql(k,j,i) = qr(k,j,i) * flag
4247                   ENDIF
4248                ELSEIF ( microphysics_morrison )  THEN
4249                   ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
4250                ELSE
4251                   IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
4252                      qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
4253                      ql(k,j,i) = qc(k,j,i) * flag
4254                   ELSE
4255                      qc(k,j,i) = 0.0_wp
4256                      ql(k,j,i) = 0.0_wp
4257                   ENDIF
4258                ENDIF
4259             ENDDO
4260          ENDDO
4261       ENDDO
4262
4263    END SUBROUTINE calc_liquid_water_content
4264
4265!------------------------------------------------------------------------------!
4266! Description:
4267! ------------
4268!> This function computes the gamma function (Press et al., 1992).
4269!> The gamma function is needed for the calculation of the evaporation
4270!> of rain drops.
4271!------------------------------------------------------------------------------!
4272    FUNCTION gamm( xx )
4273
4274       IMPLICIT NONE
4275
4276       INTEGER(iwp) ::  j            !<
4277
4278       REAL(wp)     ::  gamm         !<
4279       REAL(wp)     ::  ser          !<
4280       REAL(wp)     ::  tmp          !<
4281       REAL(wp)     ::  x_gamm       !<
4282       REAL(wp)     ::  xx           !<
4283       REAL(wp)     ::  y_gamm       !<
4284
4285
4286       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
4287       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
4288                                           -86.50532032941677_wp,      &
4289                                            24.01409824083091_wp,      &
4290                                            -1.231739572450155_wp,     &
4291                                             0.1208650973866179E-2_wp, &
4292                                            -0.5395239384953E-5_wp /)     !<
4293
4294       x_gamm = xx
4295       y_gamm = x_gamm
4296       tmp = x_gamm + 5.5_wp
4297       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
4298       ser = 1.000000000190015_wp
4299
4300       DO  j = 1, 6
4301          y_gamm = y_gamm + 1.0_wp
4302          ser    = ser + cof( j ) / y_gamm
4303       ENDDO
4304
4305!
4306!--    Until this point the algorithm computes the logarithm of the gamma
4307!--    function. Hence, the exponential function is used.
4308!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
4309       gamm = EXP( tmp ) * stp * ser / x_gamm
4310
4311       RETURN
4312
4313    END FUNCTION gamm
4314
4315 END MODULE bulk_cloud_model_mod
Note: See TracBrowser for help on using the repository browser.