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

Last change on this file since 3371 was 3371, checked in by knoop, 3 years ago

Bugfix: error in placements of CPP statement nopointer fixed

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