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

Last change on this file since 3445 was 3445, checked in by schwenkel, 3 years ago

Minor bugfix and use of subroutine for supersaturation calculation in case of cache version

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