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

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