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

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

Bugfix for profiles output

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