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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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