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

Last change on this file since 4289 was 4289, checked in by knoop, 23 months ago

Removed parameters precipitation and precipitation_amount_interval from bulk cloud model namelist

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