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

Last change on this file since 3560 was 3507, checked in by schwenkel, 5 years ago

Minor bugfixes for bcm cache version and initializing prr

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