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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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