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

Last change on this file since 3869 was 3869, checked in by knoop, 2 years ago

moving the furniture in the BCM around ;-)

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