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

Last change on this file since 3776 was 3767, checked in by raasch, 5 years ago

unused variables removed from rrd-subroutines parameter list

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