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

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

Implemented non_transport_physics module interfaces

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