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

Last change on this file since 4457 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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