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

Last change on this file since 3724 was 3724, checked in by kanani, 3 years ago

Correct double-used log_point_s units (bulk_cloud_model_mod, time_integration, turbulence_closure_mod)

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