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

Last change on this file since 4178 was 4168, checked in by suehring, 2 years ago

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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