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

Last change on this file since 3931 was 3931, checked in by schwenkel, 2 years ago

renamed module interface from non_transport_physics to non_advective_processes

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