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

Last change on this file since 3720 was 3700, checked in by knoop, 5 years ago

Moved user_define_netdf_grid into user_module.f90
Added module interface for the definition of additional timeseries

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