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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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