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

Last change on this file since 4027 was 4027, checked in by schwenkel, 5 years ago

Added microphyics scheme morrision_no_rain

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