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

Last change on this file since 4268 was 4268, checked in by schwenkel, 2 years ago

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

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