Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4495

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