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

Last change on this file since 4281 was 4268, checked in by schwenkel, 19 months ago

Introducing module interface for boundary conditions and move module specific boundary conditions into their modules

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