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

Last change on this file since 4370 was 4370, checked in by raasch, 3 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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