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

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

Moving prognostic equations of bcm into bulk_cloud_model_mod

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