Home

Context Navigation

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

Last change on this file since 3294 was 3274, checked in by knoop, 6 years ago
Modularization of all bulk cloud physics code components
Property svn:keywords set to `Id`
File size: 163.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |