Home

Context Navigation

chemistry_model_mod.f90 @ 4182

Last change on this file since 4182 was 4182, checked in by scharf, 4 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 233.0 KB

Line
1	!> @file chemistry_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 2017-2019 Leibniz Universitaet Hannover
18	! Copyright 2017-2019 Karlsruhe Institute of Technology
19	! Copyright 2017-2019 Freie Universitaet Berlin
20	!------------------------------------------------------------------------------!
21	!
22	! Current revisions:
23	! -----------------
24	!
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: chemistry_model_mod.f90 4182 2019-08-22 15:20:23Z scharf $
29	! Corrected "Former revisions" section
30	!
31	! 4167 2019-08-16 11:01:48Z suehring
32	! Changed behaviour of masked output over surface to follow terrain and ignore
33	! buildings (J.Resler, T.Gronemeier)
34	!
35	!
36	! 4166 2019-08-16 07:54:21Z resler
37	! Bugfix in decycling
38	!
39	! 4115 2019-07-24 12:50:49Z suehring
40	! Fix faulty IF statement in decycling initialization
41	!
42	! 4110 2019-07-22 17:05:21Z suehring
43	! - Decycling boundary conditions are only set at the ghost points not on the
44	! prognostic grid points
45	! - Allocation and initialization of special advection flags cs_advc_flags_s
46	! used for chemistry. These are exclusively used for chemical species in
47	! order to distinguish from the usually-used flags which might be different
48	! when decycling is applied in combination with cyclic boundary conditions.
49	! Moreover, cs_advc_flags_s considers extended zones around buildings where
50	! first-order upwind scheme is applied for the horizontal advection terms,
51	! in order to overcome high concentration peaks due to stationary numerical
52	! oscillations caused by horizontal advection discretization.
53	!
54	! 4109 2019-07-22 17:00:34Z suehring
55	! Slightly revise setting of boundary conditions at horizontal walls, use
56	! data-structure offset index instead of pre-calculate it for each facing
57	!
58	! 4080 2019-07-09 18:17:37Z suehring
59	! Restore accidantly removed limitation to positive values
60	!
61	! 4079 2019-07-09 18:04:41Z suehring
62	! Application of monotonic flux limiter for the vertical scalar advection
63	! up to the topography top (only for the cache-optimized version at the
64	! moment).
65	!
66	! 4069 2019-07-01 14:05:51Z Giersch
67	! Masked output running index mid has been introduced as a local variable to
68	! avoid runtime error (Loop variable has been modified) in time_integration
69	!
70	! 4029 2019-06-14 14:04:35Z raasch
71	! nest_chemistry option removed
72	!
73	! 4004 2019-05-24 11:32:38Z suehring
74	! in subroutine chem_parin check emiss_lod / mod_emis only
75	! when emissions_anthropogenic is activated in namelist (E.C. Chan)
76	!
77	! 3968 2019-05-13 11:04:01Z suehring
78	! - added "emiss_lod" which serves the same function as "mode_emis"
79	! both will be synchronized with emiss_lod having pirority over
80	! mode_emis (see informational messages)
81	! - modified existing error message and introduced new informational messages
82	! - CM0436 - now also applies to invalid LOD settings
83	! - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
84	! - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
85	!
86	! 3930 2019-04-24 14:57:18Z forkel
87	! Changed chem_non_transport_physics to chem_non_advective_processes
88	!
89	!
90	! 3929 2019-04-24 12:52:08Z banzhafs
91	! Correct/complete module_interface introduction for chemistry model
92	! Add subroutine chem_exchange_horiz_bounds
93	! Bug fix deposition
94	!
95	! 3784 2019-03-05 14:16:20Z banzhafs
96	! 2D output of emission fluxes
97	!
98	! 3784 2019-03-05 14:16:20Z banzhafs
99	! Bugfix, uncomment erroneous commented variable used for dry deposition.
100	! Bugfix in 3D emission output.
101	!
102	! 3784 2019-03-05 14:16:20Z banzhafs
103	! Changes related to global restructuring of location messages and introduction
104	! of additional debug messages
105	!
106	! 3784 2019-03-05 14:16:20Z banzhafs
107	! some formatting of the deposition code
108	!
109	! 3784 2019-03-05 14:16:20Z banzhafs
110	! some formatting
111	!
112	! 3784 2019-03-05 14:16:20Z banzhafs
113	! added cs_mech to USE chem_gasphase_mod
114	!
115	! 3784 2019-03-05 14:16:20Z banzhafs
116	! renamed get_mechanismname to get_mechanism_name
117	! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
118	!
119	! 3784 2019-03-05 14:16:20Z banzhafs
120	! Unused variables removed/taken care of
121	!
122	!
123	! 3784 2019-03-05 14:16:20Z forkel
124	! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
125	!
126	!
127	! 3783 2019-03-05 13:23:50Z forkel
128	! Removed forgotte write statements an some unused variables (did not touch the
129	! parts related to deposition)
130	!
131	!
132	! 3780 2019-03-05 11:19:45Z forkel
133	! Removed READ from unit 10, added CALL get_mechanismname
134	!
135	!
136	! 3767 2019-02-27 08:18:02Z raasch
137	! unused variable for file index removed from rrd-subroutines parameter list
138	!
139	! 3738 2019-02-12 17:00:45Z suehring
140	! Clean-up debug prints
141	!
142	! 3737 2019-02-12 16:57:06Z suehring
143	! Enable mesoscale offline nesting for chemistry variables as well as
144	! initialization of chemistry via dynamic input file.
145	!
146	! 3719 2019-02-06 13:10:18Z kanani
147	! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call
148	! to prognostic_equations for better overview
149	!
150	! 3700 2019-01-26 17:03:42Z knoop
151	! Some interface calls moved to module_interface + cleanup
152	!
153	! 3664 2019-01-09 14:00:37Z forkel
154	! Replaced misplaced location message by @todo
155	!
156	!
157	! 3654 2019-01-07 16:31:57Z suehring
158	! Disable misplaced location message
159	!
160	! 3652 2019-01-07 15:29:59Z forkel
161	! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
162	!
163	! 2718 2018-01-02 08:49:38Z maronga
164	! Initial revision
165	!
166	!
167	!
168	!
169	! Authors:
170	! --------
171	! @author Renate Forkel
172	! @author Farah Kanani-Suehring
173	! @author Klaus Ketelsen
174	! @author Basit Khan
175	! @author Sabine Banzhaf
176	!
177	!
178	!------------------------------------------------------------------------------!
179	! Description:
180	! ------------
181	!> Chemistry model for PALM-4U
182	!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
183	!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
184	!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
185	!> allowed to use the chemistry model in a precursor run and additionally
186	!> not using it in a main run
187	!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
188	!> @todo Separate boundary conditions for each chem spcs to be implemented
189	!> @todo Currently only total concentration are calculated. Resolved, parameterized
190	!> and chemistry fluxes although partially and some completely coded but
191	!> are not operational/activated in this version. bK.
192	!> @todo slight differences in passive scalar and chem spcs when chem reactions
193	!> turned off. Need to be fixed. bK
194	!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
195	!> @todo chemistry error messages
196	!
197	!------------------------------------------------------------------------------!
198
199	MODULE chemistry_model_mod
200
201	USE advec_s_pw_mod, &
202	ONLY: advec_s_pw
203
204	USE advec_s_up_mod, &
205	ONLY: advec_s_up
206
207	USE advec_ws, &
208	ONLY: advec_s_ws, ws_init_flags_scalar
209
210	USE diffusion_s_mod, &
211	ONLY: diffusion_s
212
213	USE kinds, &
214	ONLY: iwp, wp
215
216	USE indices, &
217	ONLY: advc_flags_s, &
218	nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt, wall_flags_0
219
220	USE pegrid, &
221	ONLY: myid, threads_per_task
222
223	USE bulk_cloud_model_mod, &
224	ONLY: bulk_cloud_model
225
226	USE control_parameters, &
227	ONLY: bc_lr_cyc, bc_ns_cyc, &
228	bc_dirichlet_l, &
229	bc_dirichlet_n, &
230	bc_dirichlet_r, &
231	bc_dirichlet_s, &
232	bc_radiation_l, &
233	bc_radiation_n, &
234	bc_radiation_r, &
235	bc_radiation_s, &
236	debug_output, &
237	dt_3d, humidity, initializing_actions, message_string, &
238	omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max, &
239	max_pr_user, &
240	monotonic_limiter_z, &
241	scalar_advec, &
242	timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
243
244	USE arrays_3d, &
245	ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
246
247	USE chem_gasphase_mod, &
248	ONLY: atol, chem_gasphase_integrate, cs_mech, get_mechanism_name, nkppctrl, &
249	nmaxfixsteps, nphot, nreact, nspec, nvar, phot_names, rtol, spc_names, t_steps, vl_dim
250
251	USE chem_modules
252
253	USE chem_photolysis_mod, &
254	ONLY: photolysis_control
255
256	USE cpulog, &
257	ONLY: cpu_log, log_point_s
258
259	USE statistics
260
261	USE surface_mod, &
262	ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
263
264	IMPLICIT NONE
265
266	PRIVATE
267	SAVE
268
269	REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_1 !< pointer for swapping of timelevels for conc
270	REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_2 !< pointer for swapping of timelevels for conc
271	REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_3 !< pointer for swapping of timelevels for conc
272	REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_av !< averaged concentrations of chemical species
273	REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: freq_1 !< pointer for phtolysis frequncies
274	!< (only 1 timelevel required)
275	INTEGER, DIMENSION(nkppctrl) :: icntrl !< 20 integer parameters for fine tuning KPP code
276	!< (e.g. solver type)
277	REAL(kind=wp), DIMENSION(nkppctrl) :: rcntrl !< 20 real parameters for fine tuning of KPP code
278	!< (e.g starting internal timestep of solver)
279	!
280	!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
281	!-- approproate for chemicals compounds since they may accumulate too much.
282	!-- If no proper boundary conditions from a DYNAMIC input file are available,
283	!-- de-cycling applies the initial profiles at the inflow boundaries for
284	!-- Dirichlet boundary conditions
285	LOGICAL :: decycle_chem_lr = .FALSE. !< switch for setting decycling in left-right direction
286	LOGICAL :: decycle_chem_ns = .FALSE. !< switch for setting decycling in south-norht direction
287	CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = &
288	(/'dirichlet','dirichlet','dirichlet','dirichlet'/)
289	!< Decycling method at horizontal boundaries,
290	!< 1=left, 2=right, 3=south, 4=north
291	!< dirichlet = initial size distribution and
292	!< chemical composition set for the ghost and
293	!< first three layers
294	!< neumann = zero gradient
295
296	REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp
297
298	!
299	!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
300	!
301	INTEGER(iwp), PARAMETER :: nlu_dep = 15 !< Number of DEPAC landuse classes (lu's)
302	INTEGER(iwp), PARAMETER :: ncmp = 10 !< Number of DEPAC gas components
303	INTEGER(iwp), PARAMETER :: nposp = 69 !< Number of possible species for deposition
304	!
305	!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1
306	INTEGER(iwp) :: ilu_grass = 1
307	INTEGER(iwp) :: ilu_arable = 2
308	INTEGER(iwp) :: ilu_permanent_crops = 3
309	INTEGER(iwp) :: ilu_coniferous_forest = 4
310	INTEGER(iwp) :: ilu_deciduous_forest = 5
311	INTEGER(iwp) :: ilu_water_sea = 6
312	INTEGER(iwp) :: ilu_urban = 7
313	INTEGER(iwp) :: ilu_other = 8
314	INTEGER(iwp) :: ilu_desert = 9
315	INTEGER(iwp) :: ilu_ice = 10
316	INTEGER(iwp) :: ilu_savanna = 11
317	INTEGER(iwp) :: ilu_tropical_forest = 12
318	INTEGER(iwp) :: ilu_water_inland = 13
319	INTEGER(iwp) :: ilu_mediterrean_scrub = 14
320	INTEGER(iwp) :: ilu_semi_natural_veg = 15
321
322	!
323	!-- NH3/SO2 ratio regimes:
324	INTEGER(iwp), PARAMETER :: iratns_low = 1 !< low ratio NH3/SO2
325	INTEGER(iwp), PARAMETER :: iratns_high = 2 !< high ratio NH3/SO2
326	INTEGER(iwp), PARAMETER :: iratns_very_low = 3 !< very low ratio NH3/SO2
327	!
328	!-- Default:
329	INTEGER, PARAMETER :: iratns_default = iratns_low
330	!
331	!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
332	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999., &
333	-999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
334	!
335	!-- Set temperatures per land use for f_temp
336	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: tmin = (/ 12.0, 12.0, 12.0, 0.0, 0.0, -999., -999., 12.0, -999., -999., &
337	12.0, 0.0, -999., 12.0, 8.0/)
338	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: topt = (/ 26.0, 26.0, 26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999., &
339	26.0, 20.0, -999., 26.0, 24.0 /)
340	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: tmax = (/ 40.0, 40.0, 40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., &
341	40.0, 35.0, -999., 40.0, 39.0 /)
342	!
343	!-- Set f_min:
344	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01, &
345	0.1, -999., 0.01, 0.04/)
346
347	!
348	!-- Set maximal conductance (m/s)
349	!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from mmol O3/m2/s to m/s
350	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999., &
351	270., 150., -999., 300., 422./)/41000.
352	!
353	!-- Set max, min for vapour pressure deficit vpd
354	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3, &
355	1.0, -999., 0.9, 2.8/)
356	REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_min = (/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999., 3.0, -999., -999., 3.0, &
357	3.25, -999., 2.8, 4.5/)
358
359	PUBLIC nreact
360	PUBLIC nspec !< number of gas phase chemical species including constant compound (e.g. N2)
361	PUBLIC nvar !< number of variable gas phase chemical species (nvar <= nspec)
362	PUBLIC spc_names !< names of gas phase chemical species (come from KPP) (come from KPP)
363	PUBLIC spec_conc_2
364	!
365	!-- Interface section
366	INTERFACE chem_3d_data_averaging
367	MODULE PROCEDURE chem_3d_data_averaging
368	END INTERFACE chem_3d_data_averaging
369
370	INTERFACE chem_boundary_conds
371	MODULE PROCEDURE chem_boundary_conds
372	MODULE PROCEDURE chem_boundary_conds_decycle
373	END INTERFACE chem_boundary_conds
374
375	INTERFACE chem_check_data_output
376	MODULE PROCEDURE chem_check_data_output
377	END INTERFACE chem_check_data_output
378
379	INTERFACE chem_data_output_2d
380	MODULE PROCEDURE chem_data_output_2d
381	END INTERFACE chem_data_output_2d
382
383	INTERFACE chem_data_output_3d
384	MODULE PROCEDURE chem_data_output_3d
385	END INTERFACE chem_data_output_3d
386
387	INTERFACE chem_data_output_mask
388	MODULE PROCEDURE chem_data_output_mask
389	END INTERFACE chem_data_output_mask
390
391	INTERFACE chem_check_data_output_pr
392	MODULE PROCEDURE chem_check_data_output_pr
393	END INTERFACE chem_check_data_output_pr
394
395	INTERFACE chem_check_parameters
396	MODULE PROCEDURE chem_check_parameters
397	END INTERFACE chem_check_parameters
398
399	INTERFACE chem_define_netcdf_grid
400	MODULE PROCEDURE chem_define_netcdf_grid
401	END INTERFACE chem_define_netcdf_grid
402
403	INTERFACE chem_header
404	MODULE PROCEDURE chem_header
405	END INTERFACE chem_header
406
407	INTERFACE chem_init_arrays
408	MODULE PROCEDURE chem_init_arrays
409	END INTERFACE chem_init_arrays
410
411	INTERFACE chem_init
412	MODULE PROCEDURE chem_init
413	END INTERFACE chem_init
414
415	INTERFACE chem_init_profiles
416	MODULE PROCEDURE chem_init_profiles
417	END INTERFACE chem_init_profiles
418
419	INTERFACE chem_integrate
420	MODULE PROCEDURE chem_integrate_ij
421	END INTERFACE chem_integrate
422
423	INTERFACE chem_parin
424	MODULE PROCEDURE chem_parin
425	END INTERFACE chem_parin
426
427	INTERFACE chem_actions
428	MODULE PROCEDURE chem_actions
429	MODULE PROCEDURE chem_actions_ij
430	END INTERFACE chem_actions
431
432	INTERFACE chem_non_advective_processes
433	MODULE PROCEDURE chem_non_advective_processes
434	MODULE PROCEDURE chem_non_advective_processes_ij
435	END INTERFACE chem_non_advective_processes
436
437	INTERFACE chem_exchange_horiz_bounds
438	MODULE PROCEDURE chem_exchange_horiz_bounds
439	END INTERFACE chem_exchange_horiz_bounds
440
441	INTERFACE chem_prognostic_equations
442	MODULE PROCEDURE chem_prognostic_equations
443	MODULE PROCEDURE chem_prognostic_equations_ij
444	END INTERFACE chem_prognostic_equations
445
446	INTERFACE chem_rrd_local
447	MODULE PROCEDURE chem_rrd_local
448	END INTERFACE chem_rrd_local
449
450	INTERFACE chem_statistics
451	MODULE PROCEDURE chem_statistics
452	END INTERFACE chem_statistics
453
454	INTERFACE chem_swap_timelevel
455	MODULE PROCEDURE chem_swap_timelevel
456	END INTERFACE chem_swap_timelevel
457
458	INTERFACE chem_wrd_local
459	MODULE PROCEDURE chem_wrd_local
460	END INTERFACE chem_wrd_local
461
462	INTERFACE chem_depo
463	MODULE PROCEDURE chem_depo
464	END INTERFACE chem_depo
465
466	INTERFACE drydepos_gas_depac
467	MODULE PROCEDURE drydepos_gas_depac
468	END INTERFACE drydepos_gas_depac
469
470	INTERFACE rc_special
471	MODULE PROCEDURE rc_special
472	END INTERFACE rc_special
473
474	INTERFACE rc_gw
475	MODULE PROCEDURE rc_gw
476	END INTERFACE rc_gw
477
478	INTERFACE rw_so2
479	MODULE PROCEDURE rw_so2
480	END INTERFACE rw_so2
481
482	INTERFACE rw_nh3_sutton
483	MODULE PROCEDURE rw_nh3_sutton
484	END INTERFACE rw_nh3_sutton
485
486	INTERFACE rw_constant
487	MODULE PROCEDURE rw_constant
488	END INTERFACE rw_constant
489
490	INTERFACE rc_gstom
491	MODULE PROCEDURE rc_gstom
492	END INTERFACE rc_gstom
493
494	INTERFACE rc_gstom_emb
495	MODULE PROCEDURE rc_gstom_emb
496	END INTERFACE rc_gstom_emb
497
498	INTERFACE par_dir_diff
499	MODULE PROCEDURE par_dir_diff
500	END INTERFACE par_dir_diff
501
502	INTERFACE rc_get_vpd
503	MODULE PROCEDURE rc_get_vpd
504	END INTERFACE rc_get_vpd
505
506	INTERFACE rc_gsoil_eff
507	MODULE PROCEDURE rc_gsoil_eff
508	END INTERFACE rc_gsoil_eff
509
510	INTERFACE rc_rinc
511	MODULE PROCEDURE rc_rinc
512	END INTERFACE rc_rinc
513
514	INTERFACE rc_rctot
515	MODULE PROCEDURE rc_rctot
516	END INTERFACE rc_rctot
517
518	! INTERFACE rc_comp_point_rc_eff
519	! MODULE PROCEDURE rc_comp_point_rc_eff
520	! END INTERFACE rc_comp_point_rc_eff
521
522	INTERFACE drydepo_aero_zhang_vd
523	MODULE PROCEDURE drydepo_aero_zhang_vd
524	END INTERFACE drydepo_aero_zhang_vd
525
526	INTERFACE get_rb_cell
527	MODULE PROCEDURE get_rb_cell
528	END INTERFACE get_rb_cell
529
530
531
532	PUBLIC chem_3d_data_averaging, chem_boundary_conds, &
533	chem_boundary_conds_decycle, chem_check_data_output, &
534	chem_check_data_output_pr, chem_check_parameters, &
535	chem_data_output_2d, chem_data_output_3d, chem_data_output_mask, &
536	chem_define_netcdf_grid, chem_header, chem_init, chem_init_arrays, &
537	chem_init_profiles, chem_integrate, chem_parin, &
538	chem_actions, chem_prognostic_equations, chem_rrd_local, &
539	chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo, &
540	chem_non_advective_processes, chem_exchange_horiz_bounds
541
542	CONTAINS
543
544
545	!------------------------------------------------------------------------------!
546	! Description:
547	! ------------
548	!> Subroutine for averaging 3D data of chemical species. Due to the fact that
549	!> the averaged chem arrays are allocated in chem_init, no if-query concerning
550	!> the allocation is required (in any mode). Attention: If you just specify an
551	!> averaged output quantity in the _p3dr file during restarts the first output
552	!> includes the time between the beginning of the restart run and the first
553	!> output time (not necessarily the whole averaging_interval you have
554	!> specified in your _p3d/_p3dr file )
555	!------------------------------------------------------------------------------!
556	SUBROUTINE chem_3d_data_averaging( mode, variable )
557
558
559	USE control_parameters
560
561	CHARACTER (LEN=*) :: mode !<
562	CHARACTER (LEN=*) :: variable !<
563
564	LOGICAL :: match_def !< flag indicating default-type surface
565	LOGICAL :: match_lsm !< flag indicating natural-type surface
566	LOGICAL :: match_usm !< flag indicating urban-type surface
567
568	INTEGER(iwp) :: i !< grid index x direction
569	INTEGER(iwp) :: j !< grid index y direction
570	INTEGER(iwp) :: k !< grid index z direction
571	INTEGER(iwp) :: m !< running index surface type
572	INTEGER(iwp) :: lsp !< running index for chem spcs
573
574	IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_') ) THEN
575
576	IF ( mode == 'allocate' ) THEN
577
578	DO lsp = 1, nspec
579	IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
580	TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN
581	chem_species(lsp)%conc_av = 0.0_wp
582	ENDIF
583	ENDDO
584
585	ELSEIF ( mode == 'sum' ) THEN
586
587	DO lsp = 1, nspec
588	IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
589	TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN
590	DO i = nxlg, nxrg
591	DO j = nysg, nyng
592	DO k = nzb, nzt+1
593	chem_species(lsp)%conc_av(k,j,i) = &
594	chem_species(lsp)%conc_av(k,j,i) + &
595	chem_species(lsp)%conc(k,j,i)
596	ENDDO
597	ENDDO
598	ENDDO
599	ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) ) THEN
600	DO i = nxl, nxr
601	DO j = nys, nyn
602	match_def = surf_def_h(0)%start_index(j,i) <= &
603	surf_def_h(0)%end_index(j,i)
604	match_lsm = surf_lsm_h%start_index(j,i) <= &
605	surf_lsm_h%end_index(j,i)
606	match_usm = surf_usm_h%start_index(j,i) <= &
607	surf_usm_h%end_index(j,i)
608
609	IF ( match_def ) THEN
610	m = surf_def_h(0)%end_index(j,i)
611	chem_species(lsp)%cssws_av(j,i) = &
612	chem_species(lsp)%cssws_av(j,i) + &
613	surf_def_h(0)%cssws(lsp,m)
614	ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN
615	m = surf_lsm_h%end_index(j,i)
616	chem_species(lsp)%cssws_av(j,i) = &
617	chem_species(lsp)%cssws_av(j,i) + &
618	surf_lsm_h%cssws(lsp,m)
619	ELSEIF ( match_usm ) THEN
620	m = surf_usm_h%end_index(j,i)
621	chem_species(lsp)%cssws_av(j,i) = &
622	chem_species(lsp)%cssws_av(j,i) + &
623	surf_usm_h%cssws(lsp,m)
624	ENDIF
625	ENDDO
626	ENDDO
627	ENDIF
628	ENDDO
629
630	ELSEIF ( mode == 'average' ) THEN
631
632	DO lsp = 1, nspec
633	IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
634	TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN
635	DO i = nxlg, nxrg
636	DO j = nysg, nyng
637	DO k = nzb, nzt+1
638	chem_species(lsp)%conc_av(k,j,i) = &
639	chem_species(lsp)%conc_av(k,j,i) / &
640	REAL( average_count_3d, KIND=wp )
641	ENDDO
642	ENDDO
643	ENDDO
644
645	ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) ) THEN
646	DO i = nxlg, nxrg
647	DO j = nysg, nyng
648	chem_species(lsp)%cssws_av(j,i) = &
649	chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
650	ENDDO
651	ENDDO
652	CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )
653	ENDIF
654	ENDDO
655	ENDIF
656
657	ENDIF
658
659	END SUBROUTINE chem_3d_data_averaging
660
661
662	!------------------------------------------------------------------------------!
663	! Description:
664	! ------------
665	!> Subroutine to initialize and set all boundary conditions for chemical species
666	!------------------------------------------------------------------------------!
667	SUBROUTINE chem_boundary_conds( mode )
668
669	USE control_parameters, &
670	ONLY: bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s
671
672	USE arrays_3d, &
673	ONLY: dzu
674
675	USE surface_mod, &
676	ONLY: bc_h
677
678	CHARACTER (LEN=*), INTENT(IN) :: mode
679	INTEGER(iwp) :: i !< grid index x direction.
680	INTEGER(iwp) :: j !< grid index y direction.
681	INTEGER(iwp) :: k !< grid index z direction.
682	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls.
683	INTEGER(iwp) :: m !< running index surface elements.
684	INTEGER(iwp) :: lsp !< running index for chem spcs.
685
686
687	SELECT CASE ( TRIM( mode ) )
688	CASE ( 'init' )
689
690	IF ( bc_cs_b == 'dirichlet' ) THEN
691	ibc_cs_b = 0
692	ELSEIF ( bc_cs_b == 'neumann' ) THEN
693	ibc_cs_b = 1
694	ELSE
695	message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'
696	CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
697	ENDIF
698	!
699	!-- Set Integer flags and check for possible erroneous settings for top
700	!-- boundary condition.
701	IF ( bc_cs_t == 'dirichlet' ) THEN
702	ibc_cs_t = 0
703	ELSEIF ( bc_cs_t == 'neumann' ) THEN
704	ibc_cs_t = 1
705	ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN
706	ibc_cs_t = 2
707	ELSEIF ( bc_cs_t == 'nested' ) THEN
708	ibc_cs_t = 3
709	ELSE
710	message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
711	CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
712	ENDIF
713
714	CASE ( 'set_bc_bottomtop' )
715	!
716	!-- Boundary condtions for chemical species at horizontal walls
717	DO lsp = 1, nspec
718	IF ( ibc_cs_b == 0 ) THEN
719	DO l = 0, 1
720	!$OMP PARALLEL DO PRIVATE( i, j, k )
721	DO m = 1, bc_h(l)%ns
722	i = bc_h(l)%i(m)
723	j = bc_h(l)%j(m)
724	k = bc_h(l)%k(m)
725	chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = &
726	chem_species(lsp)%conc(k+bc_h(l)%koff,j,i)
727	ENDDO
728	ENDDO
729
730	ELSEIF ( ibc_cs_b == 1 ) THEN
731	!
732	!-- in boundary_conds there is som extra loop over m here for passive tracer
733	DO l = 0, 1
734	!$OMP PARALLEL DO PRIVATE( i, j, k )
735	DO m = 1, bc_h(l)%ns
736	i = bc_h(l)%i(m)
737	j = bc_h(l)%j(m)
738	k = bc_h(l)%k(m)
739	chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = &
740	chem_species(lsp)%conc_p(k,j,i)
741
742	ENDDO
743	ENDDO
744	ENDIF
745	ENDDO ! end lsp loop
746	!
747	!-- Top boundary conditions for chemical species - Should this not be done for all species?
748	IF ( ibc_cs_t == 0 ) THEN
749	DO lsp = 1, nspec
750	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
751	ENDDO
752	ELSEIF ( ibc_cs_t == 1 ) THEN
753	DO lsp = 1, nspec
754	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
755	ENDDO
756	ELSEIF ( ibc_cs_t == 2 ) THEN
757	DO lsp = 1, nspec
758	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
759	ENDDO
760	ENDIF
761
762	CASE ( 'set_bc_lateral' )
763	!
764	!-- Lateral boundary conditions for chem species at inflow boundary
765	!-- are automatically set when chem_species concentration is
766	!-- initialized. The initially set value at the inflow boundary is not
767	!-- touched during time integration, hence, this boundary value remains
768	!-- at a constant value, which is the concentration that flows into the
769	!-- domain.
770	!-- Lateral boundary conditions for chem species at outflow boundary
771
772	IF ( bc_radiation_s ) THEN
773	DO lsp = 1, nspec
774	chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
775	ENDDO
776	ELSEIF ( bc_radiation_n ) THEN
777	DO lsp = 1, nspec
778	chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
779	ENDDO
780	ELSEIF ( bc_radiation_l ) THEN
781	DO lsp = 1, nspec
782	chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
783	ENDDO
784	ELSEIF ( bc_radiation_r ) THEN
785	DO lsp = 1, nspec
786	chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
787	ENDDO
788	ENDIF
789
790	END SELECT
791
792	END SUBROUTINE chem_boundary_conds
793
794
795	!------------------------------------------------------------------------------!
796	! Description:
797	! ------------
798	!> Boundary conditions for prognostic variables in chemistry: decycling in the
799	!> x-direction-
800	!> Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
801	!> approproate for chemicals compounds since they may accumulate too much.
802	!> If no proper boundary conditions from a DYNAMIC input file are available,
803	!> de-cycling applies the initial profiles at the inflow boundaries for
804	!> Dirichlet boundary conditions
805	!------------------------------------------------------------------------------!
806	SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
807
808
809	INTEGER(iwp) :: boundary !<
810	INTEGER(iwp) :: ee !<
811	INTEGER(iwp) :: copied !<
812	INTEGER(iwp) :: i !<
813	INTEGER(iwp) :: j !<
814	INTEGER(iwp) :: k !<
815	INTEGER(iwp) :: ss !<
816
817	REAL(wp), DIMENSION(nzb:nzt+1) :: cs_pr_init
818	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: cs_3d
819	REAL(wp) :: flag !< flag to mask topography grid points
820
821
822	flag = 0.0_wp
823	!
824	!-- Left and right boundaries
825	IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN
826
827	DO boundary = 1, 2
828
829	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
830	!
831	!-- Initial profile is copied to ghost and first three layers
832	ss = 1
833	ee = 0
834	IF ( boundary == 1 .AND. nxl == 0 ) THEN
835	ss = nxlg
836	ee = nxl-1
837	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
838	ss = nxr+1
839	ee = nxrg
840	ENDIF
841
842	DO i = ss, ee
843	DO j = nysg, nyng
844	DO k = nzb+1, nzt
845	flag = MERGE( 1.0_wp, 0.0_wp, &
846	BTEST( wall_flags_0(k,j,i), 0 ) )
847	cs_3d(k,j,i) = cs_pr_init(k) * flag
848	ENDDO
849	ENDDO
850	ENDDO
851
852	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
853	!
854	!-- The value at the boundary is copied to the ghost layers to simulate
855	!-- an outlet with zero gradient
856	ss = 1
857	ee = 0
858	IF ( boundary == 1 .AND. nxl == 0 ) THEN
859	ss = nxlg
860	ee = nxl-1
861	copied = nxl
862	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
863	ss = nxr+1
864	ee = nxrg
865	copied = nxr
866	ENDIF
867
868	DO i = ss, ee
869	DO j = nysg, nyng
870	DO k = nzb+1, nzt
871	flag = MERGE( 1.0_wp, 0.0_wp, &
872	BTEST( wall_flags_0(k,j,i), 0 ) )
873	cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
874	ENDDO
875	ENDDO
876	ENDDO
877
878	ELSE
879	WRITE(message_string,*) &
880	'unknown decycling method: decycle_method (', &
881	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
882	CALL message( 'chem_boundary_conds_decycle', 'CM0431', &
883	1, 2, 0, 6, 0 )
884	ENDIF
885	ENDDO
886	ENDIF
887	!
888	!-- South and north boundaries
889	IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN
890
891	DO boundary = 3, 4
892
893	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
894	!
895	!-- Initial profile is copied to ghost and first three layers
896	ss = 1
897	ee = 0
898	IF ( boundary == 3 .AND. nys == 0 ) THEN
899	ss = nysg
900	ee = nys-1
901	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
902	ss = nyn+1
903	ee = nyng
904	ENDIF
905
906	DO i = nxlg, nxrg
907	DO j = ss, ee
908	DO k = nzb+1, nzt
909	flag = MERGE( 1.0_wp, 0.0_wp, &
910	BTEST( wall_flags_0(k,j,i), 0 ) )
911	cs_3d(k,j,i) = cs_pr_init(k) * flag
912	ENDDO
913	ENDDO
914	ENDDO
915
916
917	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
918	!
919	!-- The value at the boundary is copied to the ghost layers to simulate
920	!-- an outlet with zero gradient
921	ss = 1
922	ee = 0
923	IF ( boundary == 3 .AND. nys == 0 ) THEN
924	ss = nysg
925	ee = nys-1
926	copied = nys
927	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
928	ss = nyn+1
929	ee = nyng
930	copied = nyn
931	ENDIF
932
933	DO i = nxlg, nxrg
934	DO j = ss, ee
935	DO k = nzb+1, nzt
936	flag = MERGE( 1.0_wp, 0.0_wp, &
937	BTEST( wall_flags_0(k,j,i), 0 ) )
938	cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
939	ENDDO
940	ENDDO
941	ENDDO
942
943	ELSE
944	WRITE(message_string,*) &
945	'unknown decycling method: decycle_method (', &
946	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
947	CALL message( 'chem_boundary_conds_decycle', 'CM0432', &
948	1, 2, 0, 6, 0 )
949	ENDIF
950	ENDDO
951	ENDIF
952
953
954	END SUBROUTINE chem_boundary_conds_decycle
955
956
957	!------------------------------------------------------------------------------!
958	! Description:
959	! ------------
960	!> Subroutine for checking data output for chemical species
961	!------------------------------------------------------------------------------!
962	SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
963
964
965	CHARACTER (LEN=*) :: unit !<
966	CHARACTER (LEN=*) :: var !<
967
968	INTEGER(iwp) :: i
969	INTEGER(iwp) :: lsp
970	INTEGER(iwp) :: ilen
971	INTEGER(iwp) :: k
972
973	CHARACTER(LEN=16) :: spec_name
974
975	!
976	!-- Next statement is to avoid compiler warnings about unused variables
977	IF ( ( i + ilen + k ) > 0 .OR. var(1:1) == ' ' ) CONTINUE
978
979	unit = 'illegal'
980
981	spec_name = TRIM( var(4:) ) !< var 1:3 is 'kc_' or 'em_'.
982
983	IF ( TRIM( var(1:3) ) == 'em_' ) THEN
984	DO lsp=1,nspec
985	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
986	unit = 'mol m-2 s-1'
987	ENDIF
988	!
989	!-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
990	!-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
991	!-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
992	IF (spec_name(1:2) == 'PM') THEN
993	unit = 'kg m-2 s-1'
994	ENDIF
995	ENDDO
996
997	ELSE
998
999	DO lsp=1,nspec
1000	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
1001	unit = 'ppm'
1002	ENDIF
1003	!
1004	!-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
1005	!-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1006	!-- set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
1007	IF (spec_name(1:2) == 'PM') THEN
1008	unit = 'kg m-3'
1009	ENDIF
1010	ENDDO
1011
1012	DO lsp=1,nphot
1013	IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) ) THEN
1014	unit = 'sec-1'
1015	ENDIF
1016	ENDDO
1017	ENDIF
1018
1019
1020	RETURN
1021	END SUBROUTINE chem_check_data_output
1022
1023
1024	!------------------------------------------------------------------------------!
1025	! Description:
1026	! ------------
1027	!> Subroutine for checking data output of profiles for chemistry model
1028	!------------------------------------------------------------------------------!
1029	SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1030
1031	USE arrays_3d
1032
1033	USE control_parameters, &
1034	ONLY: data_output_pr, message_string
1035
1036	USE profil_parameter
1037
1038	USE statistics
1039
1040
1041	CHARACTER (LEN=*) :: unit !<
1042	CHARACTER (LEN=*) :: variable !<
1043	CHARACTER (LEN=*) :: dopr_unit
1044	CHARACTER (LEN=16) :: spec_name
1045
1046	INTEGER(iwp) :: var_count, lsp !<
1047
1048
1049	spec_name = TRIM( variable(4:) )
1050
1051	IF ( .NOT. air_chemistry ) THEN
1052	message_string = 'data_output_pr = ' // &
1053	TRIM( data_output_pr(var_count) ) // ' is not imp' // &
1054	'lemented for air_chemistry = .FALSE.'
1055	CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
1056
1057	ELSE
1058	DO lsp = 1, nspec
1059	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
1060	cs_pr_count = cs_pr_count+1
1061	cs_pr_index(cs_pr_count) = lsp
1062	dopr_index(var_count) = pr_palm+cs_pr_count
1063	dopr_unit = 'ppm'
1064	IF (spec_name(1:2) == 'PM') THEN
1065	dopr_unit = 'kg m-3'
1066	ENDIF
1067	hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
1068	unit = dopr_unit
1069	ENDIF
1070	ENDDO
1071	ENDIF
1072
1073	END SUBROUTINE chem_check_data_output_pr
1074
1075
1076	!------------------------------------------------------------------------------!
1077	! Description:
1078	! ------------
1079	!> Check parameters routine for chemistry_model_mod
1080	!------------------------------------------------------------------------------!
1081	SUBROUTINE chem_check_parameters
1082
1083
1084	LOGICAL :: found
1085	INTEGER (iwp) :: lsp_usr !< running index for user defined chem spcs
1086	INTEGER (iwp) :: lsp !< running index for chem spcs.
1087	!
1088	!-- check for chemical reactions status
1089	IF ( chem_gasphase_on ) THEN
1090	message_string = 'Chemical reactions: ON'
1091	CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
1092	ELSEIF ( .NOT. (chem_gasphase_on) ) THEN
1093	message_string = 'Chemical reactions: OFF'
1094	CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
1095	ENDIF
1096	!
1097	!-- check for chemistry time-step
1098	IF ( call_chem_at_all_substeps ) THEN
1099	message_string = 'Chemistry is calculated at all meteorology time-step'
1100	CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
1101	ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
1102	message_string = 'Sub-time-steps are skipped for chemistry time-steps'
1103	CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
1104	ENDIF
1105	!
1106	!-- check for photolysis scheme
1107	IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant') ) THEN
1108	message_string = 'Incorrect photolysis scheme selected, please check spelling'
1109	CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
1110	ENDIF
1111	!
1112	!-- check for decycling of chem species
1113	IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) ) THEN
1114	message_string = 'Incorrect boundary conditions. Only neumann or ' &
1115	// 'dirichlet &available for decycling chemical species '
1116	CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
1117	ENDIF
1118	!
1119	!-- check for chemical mechanism used
1120	CALL get_mechanism_name
1121	IF ( chem_mechanism /= TRIM( cs_mech ) ) THEN
1122	message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
1123	CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
1124	ENDIF
1125	!
1126	!-- chem_check_parameters is called before the array chem_species is allocated!
1127	!-- temporary switch of this part of the check
1128	! RETURN !bK commented
1129
1130	CALL chem_init_internal
1131	!
1132	!-- check for initial chem species input
1133	lsp_usr = 1
1134	lsp = 1
1135	DO WHILE ( cs_name (lsp_usr) /= 'novalue')
1136	found = .FALSE.
1137	DO lsp = 1, nvar
1138	IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) ) THEN
1139	found = .TRUE.
1140	EXIT
1141	ENDIF
1142	ENDDO
1143	IF ( .NOT. found ) THEN
1144	message_string = 'Unused/incorrect input for initial surface value: ' // &
1145	TRIM( cs_name(lsp_usr) )
1146	CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 )
1147	ENDIF
1148	lsp_usr = lsp_usr + 1
1149	ENDDO
1150	!
1151	!-- check for surface emission flux chem species
1152	lsp_usr = 1
1153	lsp = 1
1154	DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
1155	found = .FALSE.
1156	DO lsp = 1, nvar
1157	IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) ) THEN
1158	found = .TRUE.
1159	EXIT
1160	ENDIF
1161	ENDDO
1162	IF ( .NOT. found ) THEN
1163	message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: ' &
1164	// TRIM( surface_csflux_name(lsp_usr) )
1165	CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 )
1166	ENDIF
1167	lsp_usr = lsp_usr + 1
1168	ENDDO
1169
1170	END SUBROUTINE chem_check_parameters
1171
1172
1173	!------------------------------------------------------------------------------!
1174	! Description:
1175	! ------------
1176	!> Subroutine defining 2D output variables for chemical species
1177	!> @todo: Remove "mode" from argument list, not used.
1178	!------------------------------------------------------------------------------!
1179	SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, &
1180	two_d, nzb_do, nzt_do, fill_value )
1181
1182
1183	CHARACTER (LEN=*) :: grid !<
1184	CHARACTER (LEN=*) :: mode !<
1185	CHARACTER (LEN=*) :: variable !<
1186	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
1187	INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb)
1188	INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1)
1189	LOGICAL :: found !<
1190	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1191	REAL(wp) :: fill_value
1192	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf
1193
1194	!
1195	!-- local variables.
1196	CHARACTER(LEN=16) :: spec_name
1197	INTEGER(iwp) :: lsp
1198	INTEGER(iwp) :: i !< grid index along x-direction
1199	INTEGER(iwp) :: j !< grid index along y-direction
1200	INTEGER(iwp) :: k !< grid index along z-direction
1201	INTEGER(iwp) :: m !< running indices for surfaces
1202	INTEGER(iwp) :: char_len !< length of a character string
1203	!
1204	!-- Next statement is to avoid compiler warnings about unused variables
1205	IF ( mode(1:1) == ' ' .OR. two_d ) CONTINUE
1206
1207	found = .FALSE.
1208	char_len = LEN_TRIM( variable )
1209
1210	spec_name = TRIM( variable(4:char_len-3) )
1211	!
1212	!-- Output of emission values, i.e. surface fluxes cssws.
1213	IF ( variable(1:3) == 'em_' ) THEN
1214
1215	local_pf = 0.0_wp
1216
1217	DO lsp = 1, nvar
1218	IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN
1219	!
1220	!-- No average output for now.
1221	DO m = 1, surf_lsm_h%ns
1222	local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) = &
1223	local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) &
1224	+ surf_lsm_h%cssws(lsp,m)
1225	ENDDO
1226	DO m = 1, surf_usm_h%ns
1227	local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) = &
1228	local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) &
1229	+ surf_usm_h%cssws(lsp,m)
1230	ENDDO
1231	grid = 'zu'
1232	found = .TRUE.
1233	ENDIF
1234	ENDDO
1235
1236	ELSE
1237
1238	DO lsp=1,nspec
1239	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) .AND. &
1240	( (variable(char_len-2:) == '_xy') .OR. &
1241	(variable(char_len-2:) == '_xz') .OR. &
1242	(variable(char_len-2:) == '_yz') ) ) THEN
1243	!
1244	!-- todo: remove or replace by "CALL message" mechanism (kanani)
1245	! IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM( variable ) // &
1246	! TRIM( chem_species(lsp)%name )
1247	IF (av == 0) THEN
1248	DO i = nxl, nxr
1249	DO j = nys, nyn
1250	DO k = nzb_do, nzt_do
1251	local_pf(i,j,k) = MERGE( &
1252	chem_species(lsp)%conc(k,j,i), &
1253	REAL( fill_value, KIND = wp ), &
1254	BTEST( wall_flags_0(k,j,i), 0 ) )
1255	ENDDO
1256	ENDDO
1257	ENDDO
1258
1259	ELSE
1260	DO i = nxl, nxr
1261	DO j = nys, nyn
1262	DO k = nzb_do, nzt_do
1263	local_pf(i,j,k) = MERGE( &
1264	chem_species(lsp)%conc_av(k,j,i), &
1265	REAL( fill_value, KIND = wp ), &
1266	BTEST( wall_flags_0(k,j,i), 0 ) )
1267	ENDDO
1268	ENDDO
1269	ENDDO
1270	ENDIF
1271	grid = 'zu'
1272	found = .TRUE.
1273	ENDIF
1274	ENDDO
1275	ENDIF
1276
1277	RETURN
1278
1279	END SUBROUTINE chem_data_output_2d
1280
1281
1282	!------------------------------------------------------------------------------!
1283	! Description:
1284	! ------------
1285	!> Subroutine defining 3D output variables for chemical species
1286	!------------------------------------------------------------------------------!
1287	SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1288
1289
1290	USE surface_mod
1291
1292	CHARACTER (LEN=*) :: variable !<
1293	INTEGER(iwp) :: av !<
1294	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1295	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1296
1297	LOGICAL :: found !<
1298
1299	REAL(wp) :: fill_value !<
1300	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf
1301	!
1302	!-- local variables
1303	CHARACTER(LEN=16) :: spec_name
1304	INTEGER(iwp) :: i
1305	INTEGER(iwp) :: j
1306	INTEGER(iwp) :: k
1307	INTEGER(iwp) :: m !< running indices for surfaces
1308	INTEGER(iwp) :: l
1309	INTEGER(iwp) :: lsp !< running index for chem spcs
1310
1311
1312	found = .FALSE.
1313	IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
1314	RETURN
1315	ENDIF
1316
1317	spec_name = TRIM( variable(4:) )
1318
1319	IF ( variable(1:3) == 'em_' ) THEN
1320
1321	DO lsp = 1, nvar !!! cssws - nvar species, chem_species - nspec species !!!
1322	IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN
1323
1324	local_pf = 0.0_wp
1325	!
1326	!-- no average for now
1327	DO m = 1, surf_usm_h%ns
1328	local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
1329	local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
1330	ENDDO
1331	DO m = 1, surf_lsm_h%ns
1332	local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
1333	local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
1334	ENDDO
1335	DO l = 0, 3
1336	DO m = 1, surf_usm_v(l)%ns
1337	local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
1338	local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
1339	ENDDO
1340	DO m = 1, surf_lsm_v(l)%ns
1341	local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
1342	local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
1343	ENDDO
1344	ENDDO
1345	found = .TRUE.
1346	ENDIF
1347	ENDDO
1348	ELSE
1349	DO lsp = 1, nspec
1350	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN
1351	!
1352	!-- todo: remove or replace by "CALL message" mechanism (kanani)
1353	! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable ) // &
1354	! TRIM( chem_species(lsp)%name )
1355	IF (av == 0) THEN
1356	DO i = nxl, nxr
1357	DO j = nys, nyn
1358	DO k = nzb_do, nzt_do
1359	local_pf(i,j,k) = MERGE( &
1360	chem_species(lsp)%conc(k,j,i), &
1361	REAL( fill_value, KIND = wp ), &
1362	BTEST( wall_flags_0(k,j,i), 0 ) )
1363	ENDDO
1364	ENDDO
1365	ENDDO
1366
1367	ELSE
1368
1369	DO i = nxl, nxr
1370	DO j = nys, nyn
1371	DO k = nzb_do, nzt_do
1372	local_pf(i,j,k) = MERGE( &
1373	chem_species(lsp)%conc_av(k,j,i),&
1374	REAL( fill_value, KIND = wp ), &
1375	BTEST( wall_flags_0(k,j,i), 0 ) )
1376	ENDDO
1377	ENDDO
1378	ENDDO
1379	ENDIF
1380	found = .TRUE.
1381	ENDIF
1382	ENDDO
1383	ENDIF
1384
1385	RETURN
1386
1387	END SUBROUTINE chem_data_output_3d
1388
1389
1390	!------------------------------------------------------------------------------!
1391	! Description:
1392	! ------------
1393	!> Subroutine defining mask output variables for chemical species
1394	!------------------------------------------------------------------------------!
1395	SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
1396
1397
1398	USE control_parameters
1399
1400	CHARACTER(LEN=16) :: spec_name
1401	CHARACTER(LEN=*) :: variable !<
1402
1403	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
1404	INTEGER(iwp) :: lsp
1405	INTEGER(iwp) :: i !< grid index along x-direction
1406	INTEGER(iwp) :: j !< grid index along y-direction
1407	INTEGER(iwp) :: k !< grid index along z-direction
1408	INTEGER(iwp) :: im !< loop index for masked variables
1409	INTEGER(iwp) :: jm !< loop index for masked variables
1410	INTEGER(iwp) :: kk !< masked output index along z-direction
1411	INTEGER(iwp) :: mid !< masked output running index
1412	INTEGER(iwp) :: ktt !< k index of highest terrain surface
1413
1414	LOGICAL :: found
1415
1416	REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
1417	local_pf !<
1418
1419	REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute
1420
1421	!
1422	!-- local variables.
1423
1424	spec_name = TRIM( variable(4:) )
1425	found = .FALSE.
1426
1427	DO lsp=1,nspec
1428	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN
1429	!
1430	!-- todo: remove or replace by "CALL message" mechanism (kanani)
1431	! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable ) // &
1432	! TRIM( chem_species(lsp)%name )
1433	IF (av == 0) THEN
1434	IF ( .NOT. mask_surface(mid) ) THEN
1435
1436	DO i = 1, mask_size_l(mid,1)
1437	DO j = 1, mask_size_l(mid,2)
1438	DO k = 1, mask_size(mid,3)
1439	local_pf(i,j,k) = chem_species(lsp)%conc( &
1440	mask_k(mid,k), &
1441	mask_j(mid,j), &
1442	mask_i(mid,i) )
1443	ENDDO
1444	ENDDO
1445	ENDDO
1446
1447	ELSE
1448	!
1449	!-- Terrain-following masked output
1450	DO i = 1, mask_size_l(mid,1)
1451	DO j = 1, mask_size_l(mid,2)
1452	!-- Get k index of the highest terraing surface
1453	im = mask_i(mid,i)
1454	jm = mask_j(mid,j)
1455	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
1456	DO k = 1, mask_size_l(mid,3)
1457	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1458	!-- Set value if not in building
1459	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
1460	local_pf(i,j,k) = fill_value
1461	ELSE
1462	local_pf(i,j,k) = chem_species(lsp)%conc(kk,jm,im)
1463	ENDIF
1464	ENDDO
1465	ENDDO
1466	ENDDO
1467
1468	ENDIF
1469	ELSE
1470	IF ( .NOT. mask_surface(mid) ) THEN
1471
1472	DO i = 1, mask_size_l(mid,1)
1473	DO j = 1, mask_size_l(mid,2)
1474	DO k = 1, mask_size_l(mid,3)
1475	local_pf(i,j,k) = chem_species(lsp)%conc_av( &
1476	mask_k(mid,k), &
1477	mask_j(mid,j), &
1478	mask_i(mid,i) )
1479	ENDDO
1480	ENDDO
1481	ENDDO
1482
1483	ELSE
1484	!
1485	!-- Terrain-following masked output
1486	DO i = 1, mask_size_l(mid,1)
1487	DO j = 1, mask_size_l(mid,2)
1488	!-- Get k index of the highest terraing surface
1489	im = mask_i(mid,i)
1490	jm = mask_j(mid,j)
1491	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
1492	DO k = 1, mask_size_l(mid,3)
1493	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1494	!-- Set value if not in building
1495	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
1496	local_pf(i,j,k) = fill_value
1497	ELSE
1498	local_pf(i,j,k) = chem_species(lsp)%conc_av(kk,jm,im)
1499	ENDIF
1500	ENDDO
1501	ENDDO
1502	ENDDO
1503
1504	ENDIF
1505
1506	ENDIF
1507	found = .TRUE.
1508	EXIT
1509	ENDIF
1510	ENDDO
1511
1512	RETURN
1513
1514	END SUBROUTINE chem_data_output_mask
1515
1516
1517	!------------------------------------------------------------------------------!
1518	! Description:
1519	! ------------
1520	!> Subroutine defining appropriate grid for netcdf variables.
1521	!> It is called out from subroutine netcdf.
1522	!------------------------------------------------------------------------------!
1523	SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1524
1525
1526	CHARACTER (LEN=*), INTENT(IN) :: var !<
1527	LOGICAL, INTENT(OUT) :: found !<
1528	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
1529	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
1530	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
1531
1532	found = .TRUE.
1533
1534	IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' ) THEN !< always the same grid for chemistry variables
1535	grid_x = 'x'
1536	grid_y = 'y'
1537	grid_z = 'zu'
1538	ELSE
1539	found = .FALSE.
1540	grid_x = 'none'
1541	grid_y = 'none'
1542	grid_z = 'none'
1543	ENDIF
1544
1545
1546	END SUBROUTINE chem_define_netcdf_grid
1547
1548
1549	!------------------------------------------------------------------------------!
1550	! Description:
1551	! ------------
1552	!> Subroutine defining header output for chemistry model
1553	!------------------------------------------------------------------------------!
1554	SUBROUTINE chem_header( io )
1555
1556
1557	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1558	INTEGER(iwp) :: lsp !< running index for chem spcs
1559	INTEGER(iwp) :: cs_fixed
1560	CHARACTER (LEN=80) :: docsflux_chr
1561	CHARACTER (LEN=80) :: docsinit_chr
1562	!
1563	! Get name of chemical mechanism from chem_gasphase_mod
1564	CALL get_mechanism_name
1565	!
1566	!-- Write chemistry model header
1567	WRITE( io, 1 )
1568	!
1569	!-- Gasphase reaction status
1570	IF ( chem_gasphase_on ) THEN
1571	WRITE( io, 2 )
1572	ELSE
1573	WRITE( io, 3 )
1574	ENDIF
1575	!
1576	!-- Chemistry time-step
1577	WRITE ( io, 4 ) cs_time_step
1578	!
1579	!-- Emission mode info
1580	!-- At the moment the evaluation is done with both emiss_lod and mode_emis
1581	!-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
1582	!-- conditions can be removed (ecc 20190513)
1583	IF ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') ) THEN
1584	WRITE ( io, 5 )
1585	ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') ) THEN
1586	WRITE ( io, 6 )
1587	ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') ) THEN
1588	WRITE ( io, 7 )
1589	ENDIF
1590	!
1591	!-- Photolysis scheme info
1592	IF ( photolysis_scheme == "simple" ) THEN
1593	WRITE( io, 8 )
1594	ELSEIF (photolysis_scheme == "constant" ) THEN
1595	WRITE( io, 9 )
1596	ENDIF
1597	!
1598	!-- Emission flux info
1599	lsp = 1
1600	docsflux_chr ='Chemical species for surface emission flux: '
1601	DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1602	docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
1603	IF ( LEN_TRIM( docsflux_chr ) >= 75 ) THEN
1604	WRITE ( io, 10 ) docsflux_chr
1605	docsflux_chr = ' '
1606	ENDIF
1607	lsp = lsp + 1
1608	ENDDO
1609
1610	IF ( docsflux_chr /= '' ) THEN
1611	WRITE ( io, 10 ) docsflux_chr
1612	ENDIF
1613	!
1614	!-- initializatoin of Surface and profile chemical species
1615
1616	lsp = 1
1617	docsinit_chr ='Chemical species for initial surface and profile emissions: '
1618	DO WHILE ( cs_name(lsp) /= 'novalue' )
1619	docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
1620	IF ( LEN_TRIM( docsinit_chr ) >= 75 ) THEN
1621	WRITE ( io, 11 ) docsinit_chr
1622	docsinit_chr = ' '
1623	ENDIF
1624	lsp = lsp + 1
1625	ENDDO
1626
1627	IF ( docsinit_chr /= '' ) THEN
1628	WRITE ( io, 11 ) docsinit_chr
1629	ENDIF
1630	!
1631	!-- number of variable and fix chemical species and number of reactions
1632	cs_fixed = nspec - nvar
1633
1634	WRITE ( io, * ) ' --> Chemical Mechanism : ', cs_mech
1635	WRITE ( io, * ) ' --> Chemical species, variable: ', nvar
1636	WRITE ( io, * ) ' --> Chemical species, fixed : ', cs_fixed
1637	WRITE ( io, * ) ' --> Total number of reactions : ', nreact
1638
1639
1640	1 FORMAT (//' Chemistry model information:'/ &
1641	' ----------------------------'/)
1642	2 FORMAT (' --> Chemical reactions are turned on')
1643	3 FORMAT (' --> Chemical reactions are turned off')
1644	4 FORMAT (' --> Time-step for chemical species: ',F6.2, ' s')
1645	5 FORMAT (' --> Emission mode = DEFAULT ')
1646	6 FORMAT (' --> Emission mode = PARAMETERIZED ')
1647	7 FORMAT (' --> Emission mode = PRE-PROCESSED ')
1648	8 FORMAT (' --> Photolysis scheme used = simple ')
1649	9 FORMAT (' --> Photolysis scheme used = constant ')
1650	10 FORMAT (/' ',A)
1651	11 FORMAT (/' ',A)
1652	!
1653	!
1654	END SUBROUTINE chem_header
1655
1656
1657	!------------------------------------------------------------------------------!
1658	! Description:
1659	! ------------
1660	!> Subroutine initializating chemistry_model_mod specific arrays
1661	!------------------------------------------------------------------------------!
1662	SUBROUTINE chem_init_arrays
1663	!
1664	!-- Please use this place to allocate required arrays
1665
1666	END SUBROUTINE chem_init_arrays
1667
1668
1669	!------------------------------------------------------------------------------!
1670	! Description:
1671	! ------------
1672	!> Subroutine initializating chemistry_model_mod
1673	!------------------------------------------------------------------------------!
1674	SUBROUTINE chem_init
1675
1676	USE chem_emissions_mod, &
1677	ONLY: chem_emissions_init
1678
1679	USE netcdf_data_input_mod, &
1680	ONLY: init_3d
1681
1682
1683	INTEGER(iwp) :: i !< running index x dimension
1684	INTEGER(iwp) :: j !< running index y dimension
1685	INTEGER(iwp) :: n !< running index for chemical species
1686
1687
1688	IF ( debug_output ) CALL debug_message( 'chem_init', 'start' )
1689	!
1690	!-- Next statement is to avoid compiler warning about unused variables
1691	IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub + &
1692	ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest + &
1693	ilu_urban ) == 0 ) CONTINUE
1694
1695	IF ( emissions_anthropogenic ) CALL chem_emissions_init
1696	!
1697	!-- Chemistry variables will be initialized if availabe from dynamic
1698	!-- input file. Note, it is possible to initialize only part of the chemistry
1699	!-- variables from dynamic input.
1700	IF ( INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
1701	DO n = 1, nspec
1702	IF ( init_3d%from_file_chem(n) ) THEN
1703	DO i = nxlg, nxrg
1704	DO j = nysg, nyng
1705	chem_species(n)%conc(:,j,i) = init_3d%chem_init(:,n)
1706	ENDDO
1707	ENDDO
1708	ENDIF
1709	ENDDO
1710	ENDIF
1711
1712	IF ( debug_output ) CALL debug_message( 'chem_init', 'end' )
1713
1714	END SUBROUTINE chem_init
1715
1716
1717	!------------------------------------------------------------------------------!
1718	! Description:
1719	! ------------
1720	!> Subroutine initializating chemistry_model_mod
1721	!> internal workaround for chem_species dependency in chem_check_parameters
1722	!------------------------------------------------------------------------------!
1723	SUBROUTINE chem_init_internal
1724
1725	USE pegrid
1726
1727	USE netcdf_data_input_mod, &
1728	ONLY: chem_emis, chem_emis_att, input_pids_dynamic, init_3d, &
1729	netcdf_data_input_chemistry_data
1730
1731	!
1732	!-- Local variables
1733	INTEGER(iwp) :: i !< running index for for horiz numerical grid points
1734	INTEGER(iwp) :: j !< running index for for horiz numerical grid points
1735	INTEGER(iwp) :: lsp !< running index for chem spcs
1736	INTEGER(iwp) :: lpr_lev !< running index for chem spcs profile level
1737
1738	IF ( emissions_anthropogenic ) THEN
1739	CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
1740	ENDIF
1741	!
1742	!-- Allocate memory for chemical species
1743	ALLOCATE( chem_species(nspec) )
1744	ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1745	ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1746	ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1747	ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1748	ALLOCATE( phot_frequen(nphot) )
1749	ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
1750	ALLOCATE( bc_cs_t_val(nspec) )
1751	!
1752	!-- Initialize arrays
1753	spec_conc_1 (:,:,:,:) = 0.0_wp
1754	spec_conc_2 (:,:,:,:) = 0.0_wp
1755	spec_conc_3 (:,:,:,:) = 0.0_wp
1756	spec_conc_av(:,:,:,:) = 0.0_wp
1757
1758
1759	DO lsp = 1, nspec
1760	chem_species(lsp)%name = spc_names(lsp)
1761
1762	chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp)
1763	chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp)
1764	chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp)
1765	chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp)
1766
1767	ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
1768	chem_species(lsp)%cssws_av = 0.0_wp
1769	!
1770	!-- The following block can be useful when emission module is not applied. &
1771	!-- if emission module is applied the following block will be overwritten.
1772	ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
1773	ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
1774	ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
1775	ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
1776	chem_species(lsp)%flux_s_cs = 0.0_wp
1777	chem_species(lsp)%flux_l_cs = 0.0_wp
1778	chem_species(lsp)%diss_s_cs = 0.0_wp
1779	chem_species(lsp)%diss_l_cs = 0.0_wp
1780	!
1781	!-- Allocate memory for initial concentration profiles
1782	!-- (concentration values come from namelist)
1783	!-- (@todo (FK): Because of this, chem_init is called in palm before
1784	!-- check_parameters, since conc_pr_init is used there.
1785	!-- We have to find another solution since chem_init should
1786	!-- eventually be called from init_3d_model!!)
1787	ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
1788	chem_species(lsp)%conc_pr_init(:) = 0.0_wp
1789
1790	ENDDO
1791	!
1792	!-- Set control flags for decycling only at lateral boundary cores, within the
1793	!-- inner cores the decycle flag is set to .False.. Even though it does not
1794	!-- affect the setting of chemistry boundary conditions, this flag is used to
1795	!-- set advection control flags appropriately.
1796	decycle_chem_lr = MERGE( decycle_chem_lr, .FALSE., &
1797	nxl == 0 .OR. nxr == nx )
1798	decycle_chem_ns = MERGE( decycle_chem_ns, .FALSE., &
1799	nys == 0 .OR. nyn == ny )
1800	!
1801	!-- For some passive scalars decycling may be enabled. This case, the lateral
1802	!-- boundary conditions are non-cyclic for these scalars (chemical species
1803	!-- and aerosols), while the other scalars may have
1804	!-- cyclic boundary conditions. However, large gradients near the boundaries
1805	!-- may produce stationary numerical oscillations near the lateral boundaries
1806	!-- when a higher-order scheme is applied near these boundaries.
1807	!-- To get rid-off this, set-up additional flags that control the order of the
1808	!-- scalar advection scheme near the lateral boundaries for passive scalars
1809	!-- with decycling.
1810	IF ( scalar_advec == 'ws-scheme' ) THEN
1811	ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1812	!
1813	!-- In case of decyling, set Neumann boundary conditions for wall_flags_0
1814	!-- bit 31 instead of cyclic boundary conditions.
1815	!-- Bit 31 is used to identify extended degradation zones (please see
1816	!-- following comment).
1817	!-- Note, since several also other modules like Salsa or other future
1818	!-- one may access this bit but may have other boundary conditions, the
1819	!-- original value of wall_flags_0 bit 31 must not be modified. Hence,
1820	!-- store the boundary conditions directly on cs_advc_flags_s.
1821	!-- cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and
1822	!-- bit 31 won't be used to control the numerical order.
1823	!-- Initialize with flag 31 only.
1824	cs_advc_flags_s = 0
1825	cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, &
1826	BTEST( wall_flags_0, 31 ) )
1827
1828	IF ( decycle_chem_ns ) THEN
1829	IF ( nys == 0 ) THEN
1830	DO i = 1, nbgp
1831	cs_advc_flags_s(:,nys-i,:) = MERGE( &
1832	IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
1833	IBCLR( cs_advc_flags_s(:,nys,:), 31 ), &
1834	BTEST( cs_advc_flags_s(:,nys,:), 31 ) &
1835	)
1836	ENDDO
1837	ENDIF
1838	IF ( nyn == ny ) THEN
1839	DO i = 1, nbgp
1840	cs_advc_flags_s(:,nyn+i,:) = MERGE( &
1841	IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
1842	IBCLR( cs_advc_flags_s(:,nyn,:), 31 ), &
1843	BTEST( cs_advc_flags_s(:,nyn,:), 31 ) &
1844	)
1845	ENDDO
1846	ENDIF
1847	ENDIF
1848	IF ( decycle_chem_lr ) THEN
1849	IF ( nxl == 0 ) THEN
1850	DO i = 1, nbgp
1851	cs_advc_flags_s(:,:,nxl-i) = MERGE( &
1852	IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
1853	IBCLR( cs_advc_flags_s(:,:,nxl), 31 ), &
1854	BTEST( cs_advc_flags_s(:,:,nxl), 31 ) &
1855	)
1856	ENDDO
1857	ENDIF
1858	IF ( nxr == nx ) THEN
1859	DO i = 1, nbgp
1860	cs_advc_flags_s(:,:,nxr+i) = MERGE( &
1861	IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
1862	IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
1863	BTEST( cs_advc_flags_s(:,:,nxr), 31 ) &
1864	)
1865	ENDDO
1866	ENDIF
1867	ENDIF
1868	!
1869	!-- To initialize advection flags appropriately, pass the boundary flags.
1870	!-- The last argument indicates that a passive scalar is treated, where
1871	!-- the horizontal advection terms are degraded already 2 grid points before
1872	!-- the lateral boundary to avoid stationary oscillations at large-gradients.
1873	!-- Also, extended degradation zones are applied, where horizontal advection of
1874	!-- passive scalars is discretized by first-order scheme at all grid points
1875	!-- that in the vicinity of buildings (<= 3 grid points). Even though no
1876	!-- building is within the numerical stencil, first-order scheme is used.
1877	!-- At fourth and fifth grid point the order of the horizontal advection scheme
1878	!-- is successively upgraded.
1879	!-- These extended degradation zones are used to avoid stationary numerical
1880	!-- oscillations, which are responsible for high concentration maxima that may
1881	!-- appear under shear-free stable conditions.
1882	CALL ws_init_flags_scalar( &
1883	bc_dirichlet_l .OR. bc_radiation_l .OR. decycle_chem_lr, &
1884	bc_dirichlet_n .OR. bc_radiation_n .OR. decycle_chem_ns, &
1885	bc_dirichlet_r .OR. bc_radiation_r .OR. decycle_chem_lr, &
1886	bc_dirichlet_s .OR. bc_radiation_s .OR. decycle_chem_ns, &
1887	cs_advc_flags_s, .TRUE. )
1888	ENDIF
1889	!
1890	!-- Initial concentration of profiles is prescribed by parameters cs_profile
1891	!-- and cs_heights in the namelist &chemistry_parameters
1892
1893	CALL chem_init_profiles
1894	!
1895	!-- In case there is dynamic input file, create a list of names for chemistry
1896	!-- initial input files. Also, initialize array that indicates whether the
1897	!-- respective variable is on file or not.
1898	IF ( input_pids_dynamic ) THEN
1899	ALLOCATE( init_3d%var_names_chem(1:nspec) )
1900	ALLOCATE( init_3d%from_file_chem(1:nspec) )
1901	init_3d%from_file_chem(:) = .FALSE.
1902
1903	DO lsp = 1, nspec
1904	init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
1905	ENDDO
1906	ENDIF
1907	!
1908	!-- Initialize model variables
1909	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1910	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1911	!
1912	!-- First model run of a possible job queue.
1913	!-- Initial profiles of the variables must be computed.
1914	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
1915	!
1916	!-- Transfer initial profiles to the arrays of the 3D model
1917	DO lsp = 1, nspec
1918	DO i = nxlg, nxrg
1919	DO j = nysg, nyng
1920	DO lpr_lev = 1, nz + 1
1921	chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
1922	ENDDO
1923	ENDDO
1924	ENDDO
1925	ENDDO
1926
1927	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
1928	THEN
1929
1930	DO lsp = 1, nspec
1931	DO i = nxlg, nxrg
1932	DO j = nysg, nyng
1933	chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init
1934	ENDDO
1935	ENDDO
1936	ENDDO
1937
1938	ENDIF
1939	!
1940	!-- If required, change the surface chem spcs at the start of the 3D run
1941	IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN
1942	DO lsp = 1, nspec
1943	chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + &
1944	cs_surface_initial_change(lsp)
1945	ENDDO
1946	ENDIF
1947	!
1948	!-- Initiale old and new time levels.
1949	DO lsp = 1, nvar
1950	chem_species(lsp)%tconc_m = 0.0_wp
1951	chem_species(lsp)%conc_p = chem_species(lsp)%conc
1952	ENDDO
1953
1954	ENDIF
1955
1956	DO lsp = 1, nphot
1957	phot_frequen(lsp)%name = phot_names(lsp)
1958	!
1959	!-- todo: remove or replace by "CALL message" mechanism (kanani)
1960	!-- IF( myid == 0 ) THEN
1961	!-- WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,TRIM( phot_names(lsp) )
1962	!-- ENDIF
1963	phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp)
1964	ENDDO
1965
1966	! CALL photolysis_init ! probably also required for restart
1967
1968	RETURN
1969
1970	END SUBROUTINE chem_init_internal
1971
1972
1973	!------------------------------------------------------------------------------!
1974	! Description:
1975	! ------------
1976	!> Subroutine defining initial vertical profiles of chemical species (given by
1977	!> namelist parameters chem_profiles and chem_heights) --> which should work
1978	!> analogue to parameters u_profile, v_profile and uv_heights)
1979	!------------------------------------------------------------------------------!
1980	SUBROUTINE chem_init_profiles
1981	!
1982	!-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
1983	!< We still need to see what has to be done in case of restart run
1984
1985	USE chem_modules
1986
1987	!
1988	!-- Local variables
1989	INTEGER :: lsp !< running index for number of species in derived data type species_def
1990	INTEGER :: lsp_usr !< running index for number of species (user defined) in cs_names, cs_profiles etc
1991	INTEGER :: lpr_lev !< running index for profile level for each chem spcs.
1992	INTEGER :: npr_lev !< the next available profile lev
1993	!
1994	!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
1995	!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
1996	!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
1997	!-- "cs_heights" are prescribed, their values will!override the constant profile given by
1998	!-- "cs_surface".
1999	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
2000	lsp_usr = 1
2001	DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) !'novalue' is the default
2002	DO lsp = 1, nspec !
2003	!
2004	!-- create initial profile (conc_pr_init) for each chemical species
2005	IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN !
2006	IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
2007	!
2008	!-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
2009	DO lpr_lev = 0, nzt+1
2010	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
2011	ENDDO
2012	ELSE
2013	IF ( cs_heights(1,1) /= 0.0_wp ) THEN
2014	message_string = 'The surface value of cs_heights must be 0.0'
2015	CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
2016	ENDIF
2017
2018	use_prescribed_profile_data = .TRUE.
2019
2020	npr_lev = 1
2021	! chem_species(lsp)%conc_pr_init(0) = 0.0_wp
2022	DO lpr_lev = 1, nz+1
2023	IF ( npr_lev < 100 ) THEN
2024	DO WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
2025	npr_lev = npr_lev + 1
2026	IF ( npr_lev == 100 ) THEN
2027	message_string = 'number of chem spcs exceeding the limit'
2028	CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )
2029	EXIT
2030	ENDIF
2031	ENDDO
2032	ENDIF
2033	IF ( npr_lev < 100 .AND. cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp ) THEN
2034	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) + &
2035	( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) / &
2036	( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) * &
2037	( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
2038	ELSE
2039	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
2040	ENDIF
2041	ENDDO
2042	ENDIF
2043	!
2044	!-- If a profile is prescribed explicity using cs_profiles and cs_heights, then
2045	!-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
2046	!-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:).
2047	ENDIF
2048	ENDDO
2049	lsp_usr = lsp_usr + 1
2050	ENDDO
2051	ENDIF
2052
2053	END SUBROUTINE chem_init_profiles
2054
2055
2056	!------------------------------------------------------------------------------!
2057	! Description:
2058	! ------------
2059	!> Subroutine to integrate chemical species in the given chemical mechanism
2060	!------------------------------------------------------------------------------!
2061	SUBROUTINE chem_integrate_ij( i, j )
2062
2063	USE statistics, &
2064	ONLY: weight_pres
2065
2066	USE control_parameters, &
2067	ONLY: dt_3d, intermediate_timestep_count, time_since_reference_point
2068
2069
2070	INTEGER,INTENT(IN) :: i
2071	INTEGER,INTENT(IN) :: j
2072	!
2073	!-- local variables
2074	INTEGER(iwp) :: lsp !< running index for chem spcs.
2075	INTEGER(iwp) :: lph !< running index for photolysis frequencies
2076	INTEGER, DIMENSION(20) :: istatus
2077	REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec) :: tmp_conc
2078	REAL(kind=wp), DIMENSION(nzb+1:nzt) :: tmp_temp
2079	REAL(kind=wp), DIMENSION(nzb+1:nzt) :: tmp_qvap
2080	REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot) :: tmp_phot
2081	REAL(kind=wp), DIMENSION(nzb+1:nzt) :: tmp_fact
2082	REAL(kind=wp), DIMENSION(nzb+1:nzt) :: tmp_fact_i !< conversion factor between
2083	!< molecules cm^{-3} and ppm
2084
2085	INTEGER,DIMENSION(nzb+1:nzt) :: nacc !< Number of accepted steps
2086	INTEGER,DIMENSION(nzb+1:nzt) :: nrej !< Number of rejected steps
2087
2088	REAL(wp) :: conv !< conversion factor
2089	REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction
2090	REAL(wp), PARAMETER :: fr2ppm = 1.0e6_wp !< Conversion factor fraction to ppm
2091	! REAL(wp), PARAMETER :: xm_air = 28.96_wp !< Mole mass of dry air
2092	! REAL(wp), PARAMETER :: xm_h2o = 18.01528_wp !< Mole mass of water vapor
2093	REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa)
2094	REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa)
2095	REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^3
2096	REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol)
2097
2098	REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
2099
2100	REAL(kind=wp) :: dt_chem
2101	!
2102	!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
2103	IF (chem_gasphase_on) THEN
2104	nacc = 0
2105	nrej = 0
2106
2107	tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
2108	!
2109	!-- convert ppm to molecules/cm**3
2110	!-- tmp_fact = 1.e-6_wp6.022e23_wp/(22.414_wp1000._wp) * 273.15_wp *
2111	!-- hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
2112	conv = ppm2fr * xna / vmolcm
2113	tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
2114	tmp_fact_i = 1.0_wp/tmp_fact
2115
2116	IF ( humidity ) THEN
2117	IF ( bulk_cloud_model ) THEN
2118	tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * &
2119	xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2120	ELSE
2121	tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2122	ENDIF
2123	ELSE
2124	tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:) !< Constant value for q if water vapor is not computed
2125	ENDIF
2126
2127	DO lsp = 1,nspec
2128	tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
2129	ENDDO
2130
2131	DO lph = 1,nphot
2132	tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
2133	ENDDO
2134	!
2135	!-- Compute length of time step
2136	IF ( call_chem_at_all_substeps ) THEN
2137	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
2138	ELSE
2139	dt_chem = dt_3d
2140	ENDIF
2141
2142	cs_time_step = dt_chem
2143
2144	IF(maxval(rcntrl) > 0.0) THEN ! Only if rcntrl is set
2145	IF( time_since_reference_point <= 2*dt_3d) THEN
2146	rcntrl_local = 0
2147	ELSE
2148	rcntrl_local = rcntrl
2149	ENDIF
2150	ELSE
2151	rcntrl_local = 0
2152	END IF
2153
2154	CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
2155	icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
2156
2157	DO lsp = 1,nspec
2158	chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
2159	ENDDO
2160
2161
2162	ENDIF
2163
2164	RETURN
2165	END SUBROUTINE chem_integrate_ij
2166
2167
2168	!------------------------------------------------------------------------------!
2169	! Description:
2170	! ------------
2171	!> Subroutine defining parin for &chemistry_parameters for chemistry model
2172	!------------------------------------------------------------------------------!
2173	SUBROUTINE chem_parin
2174
2175	USE chem_modules
2176	USE control_parameters
2177
2178	USE pegrid
2179	USE statistics
2180
2181
2182	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
2183
2184	REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping
2185	INTEGER(iwp) :: i !<
2186	INTEGER(iwp) :: max_pr_cs_tmp !<
2187
2188
2189	NAMELIST /chemistry_parameters/ bc_cs_b, &
2190	bc_cs_t, &
2191	call_chem_at_all_substeps, &
2192	chem_debug0, &
2193	chem_debug1, &
2194	chem_debug2, &
2195	chem_gasphase_on, &
2196	chem_mechanism, &
2197	cs_heights, &
2198	cs_name, &
2199	cs_profile, &
2200	cs_surface, &
2201	cs_surface_initial_change, &
2202	cs_vertical_gradient_level, &
2203	daytype_mdh, &
2204	decycle_chem_lr, &
2205	decycle_chem_ns, &
2206	decycle_method, &
2207	deposition_dry, &
2208	emissions_anthropogenic, &
2209	emiss_lod, &
2210	emiss_factor_main, &
2211	emiss_factor_side, &
2212	icntrl, &
2213	main_street_id, &
2214	max_street_id, &
2215	mode_emis, &
2216	my_steps, &
2217	rcntrl, &
2218	side_street_id, &
2219	photolysis_scheme, &
2220	wall_csflux, &
2221	cs_vertical_gradient, &
2222	top_csflux, &
2223	surface_csflux, &
2224	surface_csflux_name, &
2225	time_fac_type
2226	!
2227	!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
2228	!-- so this way we could prescribe a specific flux value for each species
2229	!> chemistry_parameters for initial profiles
2230	!> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles)
2231	!> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
2232	!> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
2233	!> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc.
2234	!> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
2235	!> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
2236	!
2237	!-- Read chem namelist
2238
2239	CHARACTER(LEN=8) :: solver_type
2240
2241	icntrl = 0
2242	rcntrl = 0.0_wp
2243	my_steps = 0.0_wp
2244	photolysis_scheme = 'simple'
2245	atol = 1.0_wp
2246	rtol = 0.01_wp
2247	!
2248	!-- Try to find chemistry package
2249	REWIND ( 11 )
2250	line = ' '
2251	DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
2252	READ ( 11, '(A)', END=20 ) line
2253	ENDDO
2254	BACKSPACE ( 11 )
2255	!
2256	!-- Read chemistry namelist
2257	READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
2258	!
2259	!-- Enable chemistry model
2260	air_chemistry = .TRUE.
2261	GOTO 20
2262
2263	10 BACKSPACE( 11 )
2264	READ( 11 , '(A)') line
2265	CALL parin_fail_message( 'chemistry_parameters', line )
2266
2267	20 CONTINUE
2268
2269	!
2270	!-- synchronize emiss_lod and mod_emis only if emissions_anthropogenic
2271	!-- is activated in the namelist. Otherwise their values are "don't care"
2272	IF ( emissions_anthropogenic ) THEN
2273	!
2274	!-- check for emission mode for chem species
2275
2276	IF ( emiss_lod < 0 ) THEN !- if LOD not defined in namelist
2277	IF ( ( mode_emis /= 'PARAMETERIZED' ) .AND. &
2278	( mode_emis /= 'DEFAULT' ) .AND. &
2279	( mode_emis /= 'PRE-PROCESSED' ) ) THEN
2280	message_string = 'Incorrect mode_emiss option select. Please check spelling'
2281	CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2282	ENDIF
2283	ELSE
2284	IF ( ( emiss_lod /= 0 ) .AND. &
2285	( emiss_lod /= 1 ) .AND. &
2286	( emiss_lod /= 2 ) ) THEN
2287	message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
2288	CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2289	ENDIF
2290	ENDIF
2291
2292	!
2293	! for reference (ecc)
2294	! IF ( (mode_emis /= 'PARAMETERIZED') .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED' ) ) THEN
2295	! message_string = 'Incorrect mode_emiss option select. Please check spelling'
2296	! CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2297	! ENDIF
2298
2299	!
2300	!-- conflict resolution for emiss_lod and mode_emis
2301	!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
2302	!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
2303	!-- this check is in place to retain backward compatibility with salsa until the
2304	!-- code is migrated completed to emiss_lod
2305	!-- note that
2306
2307	IF ( emiss_lod >= 0 ) THEN
2308
2309	SELECT CASE ( emiss_lod )
2310	CASE (0) !- parameterized mode
2311	mode_emis = 'PARAMETERIZED'
2312	CASE (1) !- default mode
2313	mode_emis = 'DEFAULT'
2314	CASE (2) !- preprocessed mode
2315	mode_emis = 'PRE-PROCESSED'
2316	END SELECT
2317
2318	message_string = 'Synchronizing mode_emis to defined emiss_lod' // &
2319	CHAR(10) // ' ' // &
2320	'NOTE - mode_emis will be depreciated in future releases' // &
2321	CHAR(10) // ' ' // &
2322	'please use emiss_lod to define emission mode'
2323
2324	CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
2325
2326	ELSE ! if emiss_lod is not set
2327
2328	SELECT CASE ( mode_emis )
2329	CASE ('PARAMETERIZED')
2330	emiss_lod = 0
2331	CASE ('DEFAULT')
2332	emiss_lod = 1
2333	CASE ('PRE-PROCESSED')
2334	emiss_lod = 2
2335	END SELECT
2336
2337	message_string = 'emiss_lod undefined. Using existing mod_emiss setting' // &
2338	CHAR(10) // ' ' // &
2339	'NOTE - mode_emis will be depreciated in future releases' // &
2340	CHAR(10) // ' ' // &
2341	' please use emiss_lod to define emission mode'
2342
2343	CALL message ( 'parin_chem', 'CM0464', 0, 0, 0, 6, 0 )
2344	ENDIF
2345
2346	ENDIF ! if emissions_anthropengic
2347
2348	t_steps = my_steps
2349	!
2350	!-- Determine the number of user-defined profiles and append them to the
2351	!-- standard data output (data_output_pr)
2352	max_pr_cs_tmp = 0
2353	i = 1
2354
2355	DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 )
2356	IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' ) THEN
2357	max_pr_cs_tmp = max_pr_cs_tmp+1
2358	ENDIF
2359	i = i +1
2360	ENDDO
2361
2362	IF ( max_pr_cs_tmp > 0 ) THEN
2363	cs_pr_namelist_found = .TRUE.
2364	max_pr_cs = max_pr_cs_tmp
2365	ENDIF
2366
2367	! Set Solver Type
2368	IF(icntrl(3) == 0) THEN
2369	solver_type = 'rodas3' !Default
2370	ELSE IF(icntrl(3) == 1) THEN
2371	solver_type = 'ros2'
2372	ELSE IF(icntrl(3) == 2) THEN
2373	solver_type = 'ros3'
2374	ELSE IF(icntrl(3) == 3) THEN
2375	solver_type = 'ro4'
2376	ELSE IF(icntrl(3) == 4) THEN
2377	solver_type = 'rodas3'
2378	ELSE IF(icntrl(3) == 5) THEN
2379	solver_type = 'rodas4'
2380	ELSE IF(icntrl(3) == 6) THEN
2381	solver_type = 'Rang3'
2382	ELSE
2383	message_string = 'illegal solver type'
2384	CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
2385	END IF
2386
2387	!
2388	!-- todo: remove or replace by "CALL message" mechanism (kanani)
2389	! write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
2390	!kk Has to be changed to right calling sequence
2391	! IF(myid == 0) THEN
2392	! write(9,*) ' '
2393	! write(9,*) 'kpp setup '
2394	! write(9,*) ' '
2395	! write(9,*) ' gas_phase chemistry: solver_type = ',TRIM( solver_type )
2396	! write(9,*) ' '
2397	! write(9,*) ' Hstart = ',rcntrl(3)
2398	! write(9,*) ' FacMin = ',rcntrl(4)
2399	! write(9,*) ' FacMax = ',rcntrl(5)
2400	! write(9,*) ' '
2401	! IF(vl_dim > 1) THEN
2402	! write(9,*) ' Vector mode vektor length = ',vl_dim
2403	! ELSE
2404	! write(9,*) ' Scalar mode'
2405	! ENDIF
2406	! write(9,*) ' '
2407	! END IF
2408
2409	RETURN
2410
2411	END SUBROUTINE chem_parin
2412
2413
2414	!------------------------------------------------------------------------------!
2415	! Description:
2416	! ------------
2417	!> Call for all grid points
2418	!------------------------------------------------------------------------------!
2419	SUBROUTINE chem_actions( location )
2420
2421
2422	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
2423
2424	SELECT CASE ( location )
2425
2426	CASE ( 'before_prognostic_equations' )
2427	!
2428	!-- Chemical reactions and deposition
2429	IF ( chem_gasphase_on ) THEN
2430	!
2431	!-- If required, calculate photolysis frequencies -
2432	!-- UNFINISHED: Why not before the intermediate timestep loop?
2433	IF ( intermediate_timestep_count == 1 ) THEN
2434	CALL photolysis_control
2435	ENDIF
2436
2437	ENDIF
2438
2439	CASE DEFAULT
2440	CONTINUE
2441
2442	END SELECT
2443
2444	END SUBROUTINE chem_actions
2445
2446
2447	!------------------------------------------------------------------------------!
2448	! Description:
2449	! ------------
2450	!> Call for grid points i,j
2451	!------------------------------------------------------------------------------!
2452
2453	SUBROUTINE chem_actions_ij( i, j, location )
2454
2455
2456	INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction
2457	INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction
2458	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
2459	INTEGER(iwp) :: dummy !< call location string
2460
2461	IF ( air_chemistry ) dummy = i + j
2462
2463	SELECT CASE ( location )
2464
2465	CASE DEFAULT
2466	CONTINUE
2467
2468	END SELECT
2469
2470
2471	END SUBROUTINE chem_actions_ij
2472
2473
2474	!------------------------------------------------------------------------------!
2475	! Description:
2476	! ------------
2477	!> Call for all grid points
2478	!------------------------------------------------------------------------------!
2479	SUBROUTINE chem_non_advective_processes()
2480
2481
2482	INTEGER(iwp) :: i !<
2483	INTEGER(iwp) :: j !<
2484
2485	!
2486	!-- Calculation of chemical reactions and deposition.
2487
2488
2489	IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps ) THEN
2490
2491	IF ( chem_gasphase_on ) THEN
2492	CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2493	!$OMP PARALLEL PRIVATE (i,j)
2494	!$OMP DO schedule(static,1)
2495	DO i = nxl, nxr
2496	DO j = nys, nyn
2497	CALL chem_integrate( i, j )
2498	ENDDO
2499	ENDDO
2500	!$OMP END PARALLEL
2501	CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2502	ENDIF
2503
2504	IF ( deposition_dry ) THEN
2505	CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2506	DO i = nxl, nxr
2507	DO j = nys, nyn
2508	CALL chem_depo( i, j )
2509	ENDDO
2510	ENDDO
2511	CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2512	ENDIF
2513
2514	ENDIF
2515
2516
2517
2518	END SUBROUTINE chem_non_advective_processes
2519
2520
2521	!------------------------------------------------------------------------------!
2522	! Description:
2523	! ------------
2524	!> Call for grid points i,j
2525	!------------------------------------------------------------------------------!
2526	SUBROUTINE chem_non_advective_processes_ij( i, j )
2527
2528
2529	INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction
2530	INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction
2531
2532	!
2533	!-- Calculation of chemical reactions and deposition.
2534
2535
2536	IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps ) THEN
2537
2538	IF ( chem_gasphase_on ) THEN
2539	CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2540	CALL chem_integrate( i, j )
2541	CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2542	ENDIF
2543
2544	IF ( deposition_dry ) THEN
2545	CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2546	CALL chem_depo( i, j )
2547	CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2548	ENDIF
2549
2550	ENDIF
2551
2552
2553
2554	END SUBROUTINE chem_non_advective_processes_ij
2555
2556	!------------------------------------------------------------------------------!
2557	! Description:
2558	! ------------
2559	!> routine for exchange horiz of chemical quantities
2560	!------------------------------------------------------------------------------!
2561	SUBROUTINE chem_exchange_horiz_bounds
2562
2563	INTEGER(iwp) :: lsp !<
2564	INTEGER(iwp) :: lsp_usr !<
2565
2566	!
2567	!-- Loop over chemical species
2568	CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
2569	DO lsp = 1, nvar
2570	CALL exchange_horiz( chem_species(lsp)%conc, nbgp )
2571	lsp_usr = 1
2572	DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
2573	IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) ) THEN
2574	!
2575	!-- As chem_exchange_horiz_bounds is called at the beginning
2576	!-- of prognostic_equations, boundary conditions are set on
2577	!-- %conc.
2578	CALL chem_boundary_conds( chem_species(lsp)%conc, &
2579	chem_species(lsp)%conc_pr_init )
2580
2581	ENDIF
2582	lsp_usr = lsp_usr + 1
2583	ENDDO
2584
2585
2586	ENDDO
2587	CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
2588
2589
2590	END SUBROUTINE chem_exchange_horiz_bounds
2591
2592
2593	!------------------------------------------------------------------------------!
2594	! Description:
2595	! ------------
2596	!> Subroutine calculating prognostic equations for chemical species
2597	!> (vector-optimized).
2598	!> Routine is called separately for each chemical species over a loop from
2599	!> prognostic_equations.
2600	!------------------------------------------------------------------------------!
2601	SUBROUTINE chem_prognostic_equations()
2602
2603
2604	INTEGER :: i !< running index
2605	INTEGER :: j !< running index
2606	INTEGER :: k !< running index
2607
2608	INTEGER(iwp) :: ilsp !<
2609
2610
2611	CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2612
2613	DO ilsp = 1, nvar
2614	!
2615	!-- Tendency terms for chemical species
2616	tend = 0.0_wp
2617	!
2618	!-- Advection terms
2619	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2620	IF ( ws_scheme_sca ) THEN
2621	CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc', &
2622	bc_dirichlet_l .OR. bc_radiation_l .OR. decycle_chem_lr, &
2623	bc_dirichlet_n .OR. bc_radiation_n .OR. decycle_chem_ns, &
2624	bc_dirichlet_r .OR. bc_radiation_r .OR. decycle_chem_lr, &
2625	bc_dirichlet_s .OR. bc_radiation_s .OR. decycle_chem_ns )
2626	ELSE
2627	CALL advec_s_pw( chem_species(ilsp)%conc )
2628	ENDIF
2629	ELSE
2630	CALL advec_s_up( chem_species(ilsp)%conc )
2631	ENDIF
2632	!
2633	!-- Diffusion terms (the last three arguments are zero)
2634	CALL diffusion_s( chem_species(ilsp)%conc, &
2635	surf_def_h(0)%cssws(ilsp,:), &
2636	surf_def_h(1)%cssws(ilsp,:), &
2637	surf_def_h(2)%cssws(ilsp,:), &
2638	surf_lsm_h%cssws(ilsp,:), &
2639	surf_usm_h%cssws(ilsp,:), &
2640	surf_def_v(0)%cssws(ilsp,:), &
2641	surf_def_v(1)%cssws(ilsp,:), &
2642	surf_def_v(2)%cssws(ilsp,:), &
2643	surf_def_v(3)%cssws(ilsp,:), &
2644	surf_lsm_v(0)%cssws(ilsp,:), &
2645	surf_lsm_v(1)%cssws(ilsp,:), &
2646	surf_lsm_v(2)%cssws(ilsp,:), &
2647	surf_lsm_v(3)%cssws(ilsp,:), &
2648	surf_usm_v(0)%cssws(ilsp,:), &
2649	surf_usm_v(1)%cssws(ilsp,:), &
2650	surf_usm_v(2)%cssws(ilsp,:), &
2651	surf_usm_v(3)%cssws(ilsp,:) )
2652	!
2653	!-- Prognostic equation for chemical species
2654	DO i = nxl, nxr
2655	DO j = nys, nyn
2656	DO k = nzb+1, nzt
2657	chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) &
2658	+ ( dt_3d * &
2659	( tsc(2) * tend(k,j,i) &
2660	+ tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) &
2661	) &
2662	- tsc(5) * rdf_sc(k) &
2663	* ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) ) &
2664	) &
2665	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2666
2667	IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp ) THEN
2668	chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)
2669	ENDIF
2670	ENDDO
2671	ENDDO
2672	ENDDO
2673	!
2674	!-- Calculate tendencies for the next Runge-Kutta step
2675	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2676	IF ( intermediate_timestep_count == 1 ) THEN
2677	DO i = nxl, nxr
2678	DO j = nys, nyn
2679	DO k = nzb+1, nzt
2680	chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2681	ENDDO
2682	ENDDO
2683	ENDDO
2684	ELSEIF ( intermediate_timestep_count < &
2685	intermediate_timestep_count_max ) THEN
2686	DO i = nxl, nxr
2687	DO j = nys, nyn
2688	DO k = nzb+1, nzt
2689	chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i) &
2690	+ 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2691	ENDDO
2692	ENDDO
2693	ENDDO
2694	ENDIF
2695	ENDIF
2696
2697	ENDDO
2698
2699	CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
2700
2701	END SUBROUTINE chem_prognostic_equations
2702
2703
2704	!------------------------------------------------------------------------------!
2705	! Description:
2706	! ------------
2707	!> Subroutine calculating prognostic equations for chemical species
2708	!> (cache-optimized).
2709	!> Routine is called separately for each chemical species over a loop from
2710	!> prognostic_equations.
2711	!------------------------------------------------------------------------------!
2712	SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
2713
2714
2715	INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
2716	INTEGER(iwp) :: ilsp
2717	!
2718	!-- local variables
2719
2720	INTEGER :: k
2721
2722	DO ilsp = 1, nvar
2723	!
2724	!-- Tendency-terms for chem spcs.
2725	tend(:,j,i) = 0.0_wp
2726	!
2727	!-- Advection terms
2728	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2729	IF ( ws_scheme_sca ) THEN
2730	CALL advec_s_ws( cs_advc_flags_s, &
2731	i, &
2732	j, &
2733	chem_species(ilsp)%conc, &
2734	'kc', &
2735	chem_species(ilsp)%flux_s_cs, &
2736	chem_species(ilsp)%diss_s_cs, &
2737	chem_species(ilsp)%flux_l_cs, &
2738	chem_species(ilsp)%diss_l_cs, &
2739	i_omp_start, &
2740	tn, &
2741	bc_dirichlet_l .OR. bc_radiation_l .OR. decycle_chem_lr, &
2742	bc_dirichlet_n .OR. bc_radiation_n .OR. decycle_chem_ns, &
2743	bc_dirichlet_r .OR. bc_radiation_r .OR. decycle_chem_lr, &
2744	bc_dirichlet_s .OR. bc_radiation_s .OR. decycle_chem_ns, &
2745	monotonic_limiter_z )
2746	ELSE
2747	CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
2748	ENDIF
2749	ELSE
2750	CALL advec_s_up( i, j, chem_species(ilsp)%conc )
2751	ENDIF
2752	!
2753	!-- Diffusion terms (the last three arguments are zero)
2754
2755	CALL diffusion_s( i, j, chem_species(ilsp)%conc, &
2756	surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), &
2757	surf_def_h(2)%cssws(ilsp,:), &
2758	surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), &
2759	surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), &
2760	surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), &
2761	surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), &
2762	surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), &
2763	surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), &
2764	surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2765	!
2766	!-- Prognostic equation for chem spcs
2767	DO k = nzb+1, nzt
2768	chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) + ( dt_3d * &
2769	( tsc(2) * tend(k,j,i) + &
2770	tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) ) &
2771	- tsc(5) * rdf_sc(k) &
2772	* ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) ) &
2773	) &
2774	* MERGE( 1.0_wp, 0.0_wp, &
2775	BTEST( wall_flags_0(k,j,i), 0 ) &
2776	)
2777
2778	IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp ) THEN
2779	chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i) !FKS6
2780	ENDIF
2781	ENDDO
2782	!
2783	!-- Calculate tendencies for the next Runge-Kutta step
2784	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2785	IF ( intermediate_timestep_count == 1 ) THEN
2786	DO k = nzb+1, nzt
2787	chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2788	ENDDO
2789	ELSEIF ( intermediate_timestep_count < &
2790	intermediate_timestep_count_max ) THEN
2791	DO k = nzb+1, nzt
2792	chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2793	5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2794	ENDDO
2795	ENDIF
2796	ENDIF
2797
2798	ENDDO
2799
2800	END SUBROUTINE chem_prognostic_equations_ij
2801
2802
2803	!------------------------------------------------------------------------------!
2804	! Description:
2805	! ------------
2806	!> Subroutine to read restart data of chemical species
2807	!------------------------------------------------------------------------------!
2808	SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
2809	nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, &
2810	nys_on_file, tmp_3d, found )
2811
2812	USE control_parameters
2813
2814
2815	CHARACTER (LEN=20) :: spc_name_av !<
2816
2817	INTEGER(iwp) :: lsp !<
2818	INTEGER(iwp) :: k !<
2819	INTEGER(iwp) :: nxlc !<
2820	INTEGER(iwp) :: nxlf !<
2821	INTEGER(iwp) :: nxl_on_file !<
2822	INTEGER(iwp) :: nxrc !<
2823	INTEGER(iwp) :: nxrf !<
2824	INTEGER(iwp) :: nxr_on_file !<
2825	INTEGER(iwp) :: nync !<
2826	INTEGER(iwp) :: nynf !<
2827	INTEGER(iwp) :: nyn_on_file !<
2828	INTEGER(iwp) :: nysc !<
2829	INTEGER(iwp) :: nysf !<
2830	INTEGER(iwp) :: nys_on_file !<
2831
2832	LOGICAL, INTENT(OUT) :: found
2833
2834	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 3D array to temp store data
2835
2836
2837	found = .FALSE.
2838
2839
2840	IF ( ALLOCATED(chem_species) ) THEN
2841
2842	DO lsp = 1, nspec
2843
2844	!< for time-averaged chemical conc.
2845	spc_name_av = TRIM( chem_species(lsp)%name )//'_av'
2846
2847	IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) ) &
2848	THEN
2849	!< read data into tmp_3d
2850	IF ( k == 1 ) READ ( 13 ) tmp_3d
2851	!< fill ..%conc in the restart run
2852	chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, &
2853	nxlc-nbgp:nxrc+nbgp) = &
2854	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2855	found = .TRUE.
2856	ELSEIF (restart_string(1:length) == spc_name_av ) THEN
2857	IF ( k == 1 ) READ ( 13 ) tmp_3d
2858	chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, &
2859	nxlc-nbgp:nxrc+nbgp) = &
2860	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2861	found = .TRUE.
2862	ENDIF
2863
2864	ENDDO
2865
2866	ENDIF
2867
2868
2869	END SUBROUTINE chem_rrd_local
2870
2871
2872	!-------------------------------------------------------------------------------!
2873	!> Description:
2874	!> Calculation of horizontally averaged profiles
2875	!> This routine is called for every statistic region (sr) defined by the user,
2876	!> but at least for the region "total domain" (sr=0).
2877	!> quantities.
2878	!-------------------------------------------------------------------------------!
2879	SUBROUTINE chem_statistics( mode, sr, tn )
2880
2881
2882	USE arrays_3d
2883
2884	USE statistics
2885
2886
2887	CHARACTER (LEN=*) :: mode !<
2888
2889	INTEGER(iwp) :: i !< running index on x-axis
2890	INTEGER(iwp) :: j !< running index on y-axis
2891	INTEGER(iwp) :: k !< vertical index counter
2892	INTEGER(iwp) :: sr !< statistical region
2893	INTEGER(iwp) :: tn !< thread number
2894	INTEGER(iwp) :: lpr !< running index chem spcs
2895
2896	IF ( mode == 'profiles' ) THEN
2897	!
2898	!
2899	!-- Sample on how to calculate horizontally averaged profiles of user-
2900	!-- defined quantities. Each quantity is identified by the index
2901	!-- "pr_palm+#" where "#" is an integer starting from 1. These
2902	!-- user-profile-numbers must also be assigned to the respective strings
2903	!-- given by data_output_pr_cs in routine user_check_data_output_pr.
2904	!-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
2905	!-- wpt), dim-4 = statistical region.
2906
2907	!$OMP DO
2908	DO i = nxl, nxr
2909	DO j = nys, nyn
2910	DO k = nzb, nzt+1
2911	DO lpr = 1, cs_pr_count
2912
2913	sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) + &
2914	chem_species(cs_pr_index(lpr))%conc(k,j,i) * &
2915	rmask(j,i,sr) * &
2916	MERGE( 1.0_wp, 0.0_wp, &
2917	BTEST( wall_flags_0(k,j,i), 22 ) )
2918	ENDDO
2919	ENDDO
2920	ENDDO
2921	ENDDO
2922	ELSEIF ( mode == 'time_series' ) THEN
2923	! @todo
2924	ENDIF
2925
2926	END SUBROUTINE chem_statistics
2927
2928
2929	!------------------------------------------------------------------------------!
2930	! Description:
2931	! ------------
2932	!> Subroutine for swapping of timelevels for chemical species
2933	!> called out from subroutine swap_timelevel
2934	!------------------------------------------------------------------------------!
2935
2936
2937	SUBROUTINE chem_swap_timelevel( level )
2938
2939
2940	INTEGER(iwp), INTENT(IN) :: level
2941	!
2942	!-- local variables
2943	INTEGER(iwp) :: lsp
2944
2945
2946	IF ( level == 0 ) THEN
2947	DO lsp=1, nvar
2948	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
2949	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
2950	ENDDO
2951	ELSE
2952	DO lsp=1, nvar
2953	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
2954	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
2955	ENDDO
2956	ENDIF
2957
2958	RETURN
2959	END SUBROUTINE chem_swap_timelevel
2960
2961
2962	!------------------------------------------------------------------------------!
2963	! Description:
2964	! ------------
2965	!> Subroutine to write restart data for chemistry model
2966	!------------------------------------------------------------------------------!
2967	SUBROUTINE chem_wrd_local
2968
2969
2970	INTEGER(iwp) :: lsp !< running index for chem spcs.
2971
2972	DO lsp = 1, nspec
2973	CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
2974	WRITE ( 14 ) chem_species(lsp)%conc
2975	CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
2976	WRITE ( 14 ) chem_species(lsp)%conc_av
2977	ENDDO
2978
2979	END SUBROUTINE chem_wrd_local
2980
2981
2982	!-------------------------------------------------------------------------------!
2983	! Description:
2984	! ------------
2985	!> Subroutine to calculate the deposition of gases and PMs. For now deposition
2986	!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
2987	!> considered. The deposition of particles is derived following Zhang et al.,
2988	!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
2989	!>
2990	!> @TODO: Consider deposition on vertical surfaces
2991	!> @TODO: Consider overlaying horizontal surfaces
2992	!> @TODO: Consider resolved vegetation
2993	!> @TODO: Check error messages
2994	!-------------------------------------------------------------------------------!
2995	SUBROUTINE chem_depo( i, j )
2996
2997	USE control_parameters, &
2998	ONLY: dt_3d, intermediate_timestep_count, latitude
2999
3000	USE arrays_3d, &
3001	ONLY: dzw, rho_air_zw
3002
3003	USE date_and_time_mod, &
3004	ONLY: day_of_year
3005
3006	USE surface_mod, &
3007	ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h, &
3008	surf_type, surf_usm_h
3009
3010	USE radiation_model_mod, &
3011	ONLY: cos_zenith
3012
3013
3014	INTEGER(iwp), INTENT(IN) :: i
3015	INTEGER(iwp), INTENT(IN) :: j
3016	INTEGER(iwp) :: k !< matching k to surface m at i,j
3017	INTEGER(iwp) :: lsp !< running index for chem spcs.
3018	INTEGER(iwp) :: luv_palm !< index of PALM LSM vegetation_type at current surface element
3019	INTEGER(iwp) :: lup_palm !< index of PALM LSM pavement_type at current surface element
3020	INTEGER(iwp) :: luw_palm !< index of PALM LSM water_type at current surface element
3021	INTEGER(iwp) :: luu_palm !< index of PALM USM walls/roofs at current surface element
3022	INTEGER(iwp) :: lug_palm !< index of PALM USM green walls/roofs at current surface element
3023	INTEGER(iwp) :: lud_palm !< index of PALM USM windows at current surface element
3024	INTEGER(iwp) :: luv_dep !< matching DEPAC LU to luv_palm
3025	INTEGER(iwp) :: lup_dep !< matching DEPAC LU to lup_palm
3026	INTEGER(iwp) :: luw_dep !< matching DEPAC LU to luw_palm
3027	INTEGER(iwp) :: luu_dep !< matching DEPAC LU to luu_palm
3028	INTEGER(iwp) :: lug_dep !< matching DEPAC LU to lug_palm
3029	INTEGER(iwp) :: lud_dep !< matching DEPAC LU to lud_palm
3030	INTEGER(iwp) :: m !< index for horizontal surfaces
3031
3032	INTEGER(iwp) :: pspec !< running index
3033	INTEGER(iwp) :: i_pspec !< index for matching depac gas component
3034	!
3035	!-- Vegetation !< Assign PALM classes to DEPAC land use classes
3036	INTEGER(iwp) :: ind_luv_user = 0 !< ERROR as no class given in PALM
3037	INTEGER(iwp) :: ind_luv_b_soil = 1 !< assigned to ilu_desert
3038	INTEGER(iwp) :: ind_luv_mixed_crops = 2 !< assigned to ilu_arable
3039	INTEGER(iwp) :: ind_luv_s_grass = 3 !< assigned to ilu_grass
3040	INTEGER(iwp) :: ind_luv_ev_needle_trees = 4 !< assigned to ilu_coniferous_forest
3041	INTEGER(iwp) :: ind_luv_de_needle_trees = 5 !< assigned to ilu_coniferous_forest
3042	INTEGER(iwp) :: ind_luv_ev_broad_trees = 6 !< assigned to ilu_tropical_forest
3043	INTEGER(iwp) :: ind_luv_de_broad_trees = 7 !< assigned to ilu_deciduous_forest
3044	INTEGER(iwp) :: ind_luv_t_grass = 8 !< assigned to ilu_grass
3045	INTEGER(iwp) :: ind_luv_desert = 9 !< assigned to ilu_desert
3046	INTEGER(iwp) :: ind_luv_tundra = 10 !< assigned to ilu_other
3047	INTEGER(iwp) :: ind_luv_irr_crops = 11 !< assigned to ilu_arable
3048	INTEGER(iwp) :: ind_luv_semidesert = 12 !< assigned to ilu_other
3049	INTEGER(iwp) :: ind_luv_ice = 13 !< assigned to ilu_ice
3050	INTEGER(iwp) :: ind_luv_marsh = 14 !< assigned to ilu_other
3051	INTEGER(iwp) :: ind_luv_ev_shrubs = 15 !< assigned to ilu_mediterrean_scrub
3052	INTEGER(iwp) :: ind_luv_de_shrubs = 16 !< assigned to ilu_mediterrean_scrub
3053	INTEGER(iwp) :: ind_luv_mixed_forest = 17 !< assigned to ilu_coniferous_forest (ave(decid+conif))
3054	INTEGER(iwp) :: ind_luv_intrup_forest = 18 !< assigned to ilu_other (ave(other+decid))
3055	!
3056	!-- Water
3057	INTEGER(iwp) :: ind_luw_user = 0 !< ERROR as no class given in PALM
3058	INTEGER(iwp) :: ind_luw_lake = 1 !< assigned to ilu_water_inland
3059	INTEGER(iwp) :: ind_luw_river = 2 !< assigned to ilu_water_inland
3060	INTEGER(iwp) :: ind_luw_ocean = 3 !< assigned to ilu_water_sea
3061	INTEGER(iwp) :: ind_luw_pond = 4 !< assigned to ilu_water_inland
3062	INTEGER(iwp) :: ind_luw_fountain = 5 !< assigned to ilu_water_inland
3063	!
3064	!-- Pavement
3065	INTEGER(iwp) :: ind_lup_user = 0 !< ERROR as no class given in PALM
3066	INTEGER(iwp) :: ind_lup_asph_conc = 1 !< assigned to ilu_desert
3067	INTEGER(iwp) :: ind_lup_asph = 2 !< assigned to ilu_desert
3068	INTEGER(iwp) :: ind_lup_conc = 3 !< assigned to ilu_desert
3069	INTEGER(iwp) :: ind_lup_sett = 4 !< assigned to ilu_desert
3070	INTEGER(iwp) :: ind_lup_pav_stones = 5 !< assigned to ilu_desert
3071	INTEGER(iwp) :: ind_lup_cobblest = 6 !< assigned to ilu_desert
3072	INTEGER(iwp) :: ind_lup_metal = 7 !< assigned to ilu_desert
3073	INTEGER(iwp) :: ind_lup_wood = 8 !< assigned to ilu_desert
3074	INTEGER(iwp) :: ind_lup_gravel = 9 !< assigned to ilu_desert
3075	INTEGER(iwp) :: ind_lup_f_gravel = 10 !< assigned to ilu_desert
3076	INTEGER(iwp) :: ind_lup_pebblest = 11 !< assigned to ilu_desert
3077	INTEGER(iwp) :: ind_lup_woodchips = 12 !< assigned to ilu_desert
3078	INTEGER(iwp) :: ind_lup_tartan = 13 !< assigned to ilu_desert
3079	INTEGER(iwp) :: ind_lup_art_turf = 14 !< assigned to ilu_desert
3080	INTEGER(iwp) :: ind_lup_clay = 15 !< assigned to ilu_desert
3081	!
3082	!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
3083	INTEGER(iwp) :: ind_p_size = 1 !< index for partsize in particle_pars
3084	INTEGER(iwp) :: ind_p_dens = 2 !< index for rhopart in particle_pars
3085	INTEGER(iwp) :: ind_p_slip = 3 !< index for slipcor in particle_pars
3086
3087	INTEGER(iwp) :: part_type !< index for particle type (PM10 or PM25) in particle_pars
3088
3089	INTEGER(iwp) :: nwet !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
3090
3091	REAL(wp) :: dt_chem !< length of chem time step
3092	REAL(wp) :: dh !< vertical grid size
3093	REAL(wp) :: inv_dh !< inverse of vertical grid size
3094	REAL(wp) :: dt_dh !< dt_chem/dh
3095
3096	REAL(wp) :: dens !< density at layer k at i,j
3097	REAL(wp) :: r_aero_surf !< aerodynamic resistance (s/m) at current surface element
3098	REAL(wp) :: ustar_surf !< ustar at current surface element
3099	REAL(wp) :: z0h_surf !< roughness length for heat at current surface element
3100	REAL(wp) :: solar_rad !< solar radiation, direct and diffuse, at current surface element
3101	REAL(wp) :: ppm2ugm3 !< conversion factor from ppm to ug/m3
3102	REAL(wp) :: rh_surf !< relative humidity at current surface element
3103	REAL(wp) :: lai !< leaf area index at current surface element
3104	REAL(wp) :: sai !< surface area index at current surface element assumed to be lai + 1
3105
3106	REAL(wp) :: slinnfac
3107	REAL(wp) :: visc !< Viscosity
3108	REAL(wp) :: vs !< Sedimentation velocity
3109	REAL(wp) :: vd_lu !< deposition velocity (m/s)
3110	REAL(wp) :: rs !< Sedimentaion resistance (s/m)
3111	REAL(wp) :: rb !< quasi-laminar boundary layer resistance (s/m)
3112	REAL(wp) :: rc_tot !< total canopy resistance (s/m)
3113
3114	REAL(wp) :: conc_ijk_ugm3 !< concentration at i, j, k in ug/m3
3115	REAL(wp) :: diffusivity !< diffusivity
3116
3117
3118	REAL(wp), DIMENSION(nspec) :: bud_luv !< budget for LSM vegetation type at current surface element
3119	REAL(wp), DIMENSION(nspec) :: bud_lup !< budget for LSM pavement type at current surface element
3120	REAL(wp), DIMENSION(nspec) :: bud_luw !< budget for LSM water type at current surface element
3121	REAL(wp), DIMENSION(nspec) :: bud_luu !< budget for USM walls/roofs at current surface element
3122	REAL(wp), DIMENSION(nspec) :: bud_lug !< budget for USM green surfaces at current surface element
3123	REAL(wp), DIMENSION(nspec) :: bud_lud !< budget for USM windows at current surface element
3124	REAL(wp), DIMENSION(nspec) :: bud !< overall budget at current surface element
3125	REAL(wp), DIMENSION(nspec) :: conc_ijk !< concentration at i,j,k
3126	REAL(wp), DIMENSION(nspec) :: ccomp_tot !< total compensation point (ug/m3), for now kept to zero for all species!
3127
3128
3129	REAL(wp) :: temp_tmp !< temperatur at i,j,k
3130	REAL(wp) :: ts !< surface temperatur in degrees celsius
3131	REAL(wp) :: qv_tmp !< surface mixing ratio at current surface element
3132	!
3133	!-- Particle parameters (PM10 (1), PM25 (2))
3134	!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
3135	!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
3136	REAL(wp), DIMENSION(1:3,1:2), PARAMETER :: particle_pars = RESHAPE( (/ &
3137	8.0e-6_wp, 1.14e3_wp, 1.016_wp, & !< 1
3138	0.7e-6_wp, 1.14e3_wp, 1.082_wp & !< 2
3139	/), (/ 3, 2 /) )
3140
3141	LOGICAL :: match_lsm !< flag indicating natural-type surface
3142	LOGICAL :: match_usm !< flag indicating urban-type surface
3143	!
3144	!-- List of names of possible tracers
3145	CHARACTER(LEN=*), PARAMETER :: pspecnames(nposp) = (/ &
3146	'NO2 ', & !< NO2
3147	'NO ', & !< NO
3148	'O3 ', & !< O3
3149	'CO ', & !< CO
3150	'form ', & !< FORM
3151	'ald ', & !< ALD
3152	'pan ', & !< PAN
3153	'mgly ', & !< MGLY
3154	'par ', & !< PAR
3155	'ole ', & !< OLE
3156	'eth ', & !< ETH
3157	'tol ', & !< TOL
3158	'cres ', & !< CRES
3159	'xyl ', & !< XYL
3160	'SO4a_f ', & !< SO4a_f
3161	'SO2 ', & !< SO2
3162	'HNO2 ', & !< HNO2
3163	'CH4 ', & !< CH4
3164	'NH3 ', & !< NH3
3165	'NO3 ', & !< NO3
3166	'OH ', & !< OH
3167	'HO2 ', & !< HO2
3168	'N2O5 ', & !< N2O5
3169	'SO4a_c ', & !< SO4a_c
3170	'NH4a_f ', & !< NH4a_f
3171	'NO3a_f ', & !< NO3a_f
3172	'NO3a_c ', & !< NO3a_c
3173	'C2O3 ', & !< C2O3
3174	'XO2 ', & !< XO2
3175	'XO2N ', & !< XO2N
3176	'cro ', & !< CRO
3177	'HNO3 ', & !< HNO3
3178	'H2O2 ', & !< H2O2
3179	'iso ', & !< ISO
3180	'ispd ', & !< ISPD
3181	'to2 ', & !< TO2
3182	'open ', & !< OPEN
3183	'terp ', & !< TERP
3184	'ec_f ', & !< EC_f
3185	'ec_c ', & !< EC_c
3186	'pom_f ', & !< POM_f
3187	'pom_c ', & !< POM_c
3188	'ppm_f ', & !< PPM_f
3189	'ppm_c ', & !< PPM_c
3190	'na_ff ', & !< Na_ff
3191	'na_f ', & !< Na_f
3192	'na_c ', & !< Na_c
3193	'na_cc ', & !< Na_cc
3194	'na_ccc ', & !< Na_ccc
3195	'dust_ff ', & !< dust_ff
3196	'dust_f ', & !< dust_f
3197	'dust_c ', & !< dust_c
3198	'dust_cc ', & !< dust_cc
3199	'dust_ccc ', & !< dust_ccc
3200	'tpm10 ', & !< tpm10
3201	'tpm25 ', & !< tpm25
3202	'tss ', & !< tss
3203	'tdust ', & !< tdust
3204	'tc ', & !< tc
3205	'tcg ', & !< tcg
3206	'tsoa ', & !< tsoa
3207	'tnmvoc ', & !< tnmvoc
3208	'SOxa ', & !< SOxa
3209	'NOya ', & !< NOya
3210	'NHxa ', & !< NHxa
3211	'NO2_obs ', & !< NO2_obs
3212	'tpm10_biascorr', & !< tpm10_biascorr
3213	'tpm25_biascorr', & !< tpm25_biascorr
3214	'O3_biascorr ' /) !< o3_biascorr
3215	!
3216	!-- tracer mole mass:
3217	REAL(wp), PARAMETER :: specmolm(nposp) = (/ &
3218	xm_O * 2 + xm_N, & !< NO2
3219	xm_O + xm_N, & !< NO
3220	xm_O * 3, & !< O3
3221	xm_C + xm_O, & !< CO
3222	xm_H * 2 + xm_C + xm_O, & !< FORM
3223	xm_H * 3 + xm_C * 2 + xm_O, & !< ALD
3224	xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, & !< PAN
3225	xm_H * 4 + xm_C * 3 + xm_O * 2, & !< MGLY
3226	xm_H * 3 + xm_C, & !< PAR
3227	xm_H * 3 + xm_C * 2, & !< OLE
3228	xm_H * 4 + xm_C * 2, & !< ETH
3229	xm_H * 8 + xm_C * 7, & !< TOL
3230	xm_H * 8 + xm_C * 7 + xm_O, & !< CRES
3231	xm_H * 10 + xm_C * 8, & !< XYL
3232	xm_S + xm_O * 4, & !< SO4a_f
3233	xm_S + xm_O * 2, & !< SO2
3234	xm_H + xm_O * 2 + xm_N, & !< HNO2
3235	xm_H * 4 + xm_C, & !< CH4
3236	xm_H * 3 + xm_N, & !< NH3
3237	xm_O * 3 + xm_N, & !< NO3
3238	xm_H + xm_O, & !< OH
3239	xm_H + xm_O * 2, & !< HO2
3240	xm_O * 5 + xm_N * 2, & !< N2O5
3241	xm_S + xm_O * 4, & !< SO4a_c
3242	xm_H * 4 + xm_N, & !< NH4a_f
3243	xm_O * 3 + xm_N, & !< NO3a_f
3244	xm_O * 3 + xm_N, & !< NO3a_c
3245	xm_C * 2 + xm_O * 3, & !< C2O3
3246	xm_dummy, & !< XO2
3247	xm_dummy, & !< XO2N
3248	xm_dummy, & !< CRO
3249	xm_H + xm_O * 3 + xm_N, & !< HNO3
3250	xm_H * 2 + xm_O * 2, & !< H2O2
3251	xm_H * 8 + xm_C * 5, & !< ISO
3252	xm_dummy, & !< ISPD
3253	xm_dummy, & !< TO2
3254	xm_dummy, & !< OPEN
3255	xm_H * 16 + xm_C * 10, & !< TERP
3256	xm_dummy, & !< EC_f
3257	xm_dummy, & !< EC_c
3258	xm_dummy, & !< POM_f
3259	xm_dummy, & !< POM_c
3260	xm_dummy, & !< PPM_f
3261	xm_dummy, & !< PPM_c
3262	xm_Na, & !< Na_ff
3263	xm_Na, & !< Na_f
3264	xm_Na, & !< Na_c
3265	xm_Na, & !< Na_cc
3266	xm_Na, & !< Na_ccc
3267	xm_dummy, & !< dust_ff
3268	xm_dummy, & !< dust_f
3269	xm_dummy, & !< dust_c
3270	xm_dummy, & !< dust_cc
3271	xm_dummy, & !< dust_ccc
3272	xm_dummy, & !< tpm10
3273	xm_dummy, & !< tpm25
3274	xm_dummy, & !< tss
3275	xm_dummy, & !< tdust
3276	xm_dummy, & !< tc
3277	xm_dummy, & !< tcg
3278	xm_dummy, & !< tsoa
3279	xm_dummy, & !< tnmvoc
3280	xm_dummy, & !< SOxa
3281	xm_dummy, & !< NOya
3282	xm_dummy, & !< NHxa
3283	xm_O * 2 + xm_N, & !< NO2_obs
3284	xm_dummy, & !< tpm10_biascorr
3285	xm_dummy, & !< tpm25_biascorr
3286	xm_O * 3 /) !< o3_biascorr
3287	!
3288	!-- Initialize surface element m
3289	m = 0
3290	k = 0
3291	!
3292	!-- LSM or USM surface present at i,j:
3293	!-- Default surfaces are NOT considered for deposition
3294	match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
3295	match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
3296	!
3297	!--For LSM surfaces
3298
3299	IF ( match_lsm ) THEN
3300	!
3301	!-- Get surface element information at i,j:
3302	m = surf_lsm_h%start_index(j,i)
3303	k = surf_lsm_h%k(m)
3304	!
3305	!-- Get needed variables for surface element m
3306	ustar_surf = surf_lsm_h%us(m)
3307	z0h_surf = surf_lsm_h%z0h(m)
3308	r_aero_surf = surf_lsm_h%r_a(m)
3309	solar_rad = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
3310
3311	lai = surf_lsm_h%lai(m)
3312	sai = lai + 1
3313	!
3314	!-- For small grid spacing neglect R_a
3315	IF ( dzw(k) <= 1.0 ) THEN
3316	r_aero_surf = 0.0_wp
3317	ENDIF
3318	!
3319	!-- Initialize lu's
3320	luv_palm = 0
3321	luv_dep = 0
3322	lup_palm = 0
3323	lup_dep = 0
3324	luw_palm = 0
3325	luw_dep = 0
3326	!
3327	!-- Initialize budgets
3328	bud_luv = 0.0_wp
3329	bud_lup = 0.0_wp
3330	bud_luw = 0.0_wp
3331	!
3332	!-- Get land use for i,j and assign to DEPAC lu
3333	IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 ) THEN
3334	luv_palm = surf_lsm_h%vegetation_type(m)
3335	IF ( luv_palm == ind_luv_user ) THEN
3336	message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3337	CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
3338	ELSEIF ( luv_palm == ind_luv_b_soil ) THEN
3339	luv_dep = 9
3340	ELSEIF ( luv_palm == ind_luv_mixed_crops ) THEN
3341	luv_dep = 2
3342	ELSEIF ( luv_palm == ind_luv_s_grass ) THEN
3343	luv_dep = 1
3344	ELSEIF ( luv_palm == ind_luv_ev_needle_trees ) THEN
3345	luv_dep = 4
3346	ELSEIF ( luv_palm == ind_luv_de_needle_trees ) THEN
3347	luv_dep = 4
3348	ELSEIF ( luv_palm == ind_luv_ev_broad_trees ) THEN
3349	luv_dep = 12
3350	ELSEIF ( luv_palm == ind_luv_de_broad_trees ) THEN
3351	luv_dep = 5
3352	ELSEIF ( luv_palm == ind_luv_t_grass ) THEN
3353	luv_dep = 1
3354	ELSEIF ( luv_palm == ind_luv_desert ) THEN
3355	luv_dep = 9
3356	ELSEIF ( luv_palm == ind_luv_tundra ) THEN
3357	luv_dep = 8
3358	ELSEIF ( luv_palm == ind_luv_irr_crops ) THEN
3359	luv_dep = 2
3360	ELSEIF ( luv_palm == ind_luv_semidesert ) THEN
3361	luv_dep = 8
3362	ELSEIF ( luv_palm == ind_luv_ice ) THEN
3363	luv_dep = 10
3364	ELSEIF ( luv_palm == ind_luv_marsh ) THEN
3365	luv_dep = 8
3366	ELSEIF ( luv_palm == ind_luv_ev_shrubs ) THEN
3367	luv_dep = 14
3368	ELSEIF ( luv_palm == ind_luv_de_shrubs ) THEN
3369	luv_dep = 14
3370	ELSEIF ( luv_palm == ind_luv_mixed_forest ) THEN
3371	luv_dep = 4
3372	ELSEIF ( luv_palm == ind_luv_intrup_forest ) THEN
3373	luv_dep = 8
3374	ENDIF
3375	ENDIF
3376
3377	IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 ) THEN
3378	lup_palm = surf_lsm_h%pavement_type(m)
3379	IF ( lup_palm == ind_lup_user ) THEN
3380	message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3381	CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
3382	ELSEIF ( lup_palm == ind_lup_asph_conc ) THEN
3383	lup_dep = 9
3384	ELSEIF ( lup_palm == ind_lup_asph ) THEN
3385	lup_dep = 9
3386	ELSEIF ( lup_palm == ind_lup_conc ) THEN
3387	lup_dep = 9
3388	ELSEIF ( lup_palm == ind_lup_sett ) THEN
3389	lup_dep = 9
3390	ELSEIF ( lup_palm == ind_lup_pav_stones ) THEN
3391	lup_dep = 9
3392	ELSEIF ( lup_palm == ind_lup_cobblest ) THEN
3393	lup_dep = 9
3394	ELSEIF ( lup_palm == ind_lup_metal ) THEN
3395	lup_dep = 9
3396	ELSEIF ( lup_palm == ind_lup_wood ) THEN
3397	lup_dep = 9
3398	ELSEIF ( lup_palm == ind_lup_gravel ) THEN
3399	lup_dep = 9
3400	ELSEIF ( lup_palm == ind_lup_f_gravel ) THEN
3401	lup_dep = 9
3402	ELSEIF ( lup_palm == ind_lup_pebblest ) THEN
3403	lup_dep = 9
3404	ELSEIF ( lup_palm == ind_lup_woodchips ) THEN
3405	lup_dep = 9
3406	ELSEIF ( lup_palm == ind_lup_tartan ) THEN
3407	lup_dep = 9
3408	ELSEIF ( lup_palm == ind_lup_art_turf ) THEN
3409	lup_dep = 9
3410	ELSEIF ( lup_palm == ind_lup_clay ) THEN
3411	lup_dep = 9
3412	ENDIF
3413	ENDIF
3414
3415	IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 ) THEN
3416	luw_palm = surf_lsm_h%water_type(m)
3417	IF ( luw_palm == ind_luw_user ) THEN
3418	message_string = 'No water type defined. Please define water type to enable deposition calculation'
3419	CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
3420	ELSEIF ( luw_palm == ind_luw_lake ) THEN
3421	luw_dep = 13
3422	ELSEIF ( luw_palm == ind_luw_river ) THEN
3423	luw_dep = 13
3424	ELSEIF ( luw_palm == ind_luw_ocean ) THEN
3425	luw_dep = 6
3426	ELSEIF ( luw_palm == ind_luw_pond ) THEN
3427	luw_dep = 13
3428	ELSEIF ( luw_palm == ind_luw_fountain ) THEN
3429	luw_dep = 13
3430	ENDIF
3431	ENDIF
3432	!
3433	!-- Set wetness indicator to dry or wet for lsm vegetation or pavement
3434	IF ( surf_lsm_h%c_liq(m) > 0 ) THEN
3435	nwet = 1
3436	ELSE
3437	nwet = 0
3438	ENDIF
3439	!
3440	!-- Compute length of time step
3441	IF ( call_chem_at_all_substeps ) THEN
3442	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3443	ELSE
3444	dt_chem = dt_3d
3445	ENDIF
3446
3447	dh = dzw(k)
3448	inv_dh = 1.0_wp / dh
3449	dt_dh = dt_chem/dh
3450	!
3451	!-- Concentration at i,j,k
3452	DO lsp = 1, nspec
3453	conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3454	ENDDO
3455
3456	!-- Temperature at i,j,k
3457	temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3458
3459	ts = temp_tmp - 273.15 !< in degrees celcius
3460	!
3461	!-- Viscosity of air
3462	visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
3463	!
3464	!-- Air density at k
3465	dens = rho_air_zw(k)
3466	!
3467	!-- Calculate relative humidity from specific humidity for DEPAC
3468	qv_tmp = MAX(q(k,j,i),0.0_wp)
3469	rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
3470	!
3471	!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3472	!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3473	!
3474	!-- Vegetation
3475	IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 ) THEN
3476
3477	!
3478	!-- No vegetation on bare soil, desert or ice:
3479	IF ( ( luv_palm == ind_luv_b_soil ) .OR. &
3480	( luv_palm == ind_luv_desert ) .OR. &
3481	( luv_palm == ind_luv_ice ) ) THEN
3482
3483	lai = 0.0_wp
3484	sai = 0.0_wp
3485
3486	ENDIF
3487
3488	slinnfac = 1.0_wp
3489	!
3490	!-- Get deposition velocity vd
3491	DO lsp = 1, nvar
3492	!
3493	!-- Initialize
3494	vs = 0.0_wp
3495	vd_lu = 0.0_wp
3496	rs = 0.0_wp
3497	rb = 0.0_wp
3498	rc_tot = 0.0_wp
3499	IF ( spc_names(lsp) == 'PM10' ) THEN
3500	part_type = 1
3501	!
3502	!-- Sedimentation velocity
3503	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3504	particle_pars(ind_p_size, part_type), &
3505	particle_pars(ind_p_slip, part_type), &
3506	visc)
3507
3508	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3509	vs, &
3510	particle_pars(ind_p_size, part_type), &
3511	particle_pars(ind_p_slip, part_type), &
3512	nwet, temp_tmp, dens, visc, &
3513	luv_dep, &
3514	r_aero_surf, ustar_surf )
3515
3516	bud_luv(lsp) = - conc_ijk(lsp) * &
3517	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3518
3519
3520	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
3521	part_type = 2
3522	!
3523	!-- Sedimentation velocity
3524	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3525	particle_pars(ind_p_size, part_type), &
3526	particle_pars(ind_p_slip, part_type), &
3527	visc)
3528
3529	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3530	vs, &
3531	particle_pars(ind_p_size, part_type), &
3532	particle_pars(ind_p_slip, part_type), &
3533	nwet, temp_tmp, dens, visc, &
3534	luv_dep , &
3535	r_aero_surf, ustar_surf )
3536
3537	bud_luv(lsp) = - conc_ijk(lsp) * &
3538	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3539
3540	ELSE !< GASES
3541	!
3542	!-- Read spc_name of current species for gas parameter
3543	IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN
3544	i_pspec = 0
3545	DO pspec = 1, nposp
3546	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
3547	i_pspec = pspec
3548	END IF
3549	ENDDO
3550
3551	ELSE
3552	!
3553	!-- For now species not deposited
3554	CYCLE
3555	ENDIF
3556	!
3557	!-- Factor used for conversion from ppb to ug/m3 :
3558	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3559	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3560	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
3561	!-- thus:
3562	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3563	!-- Use density at k:
3564
3565	ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3
3566	!
3567	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
3568	! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
3569	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
3570	!
3571	!-- Diffusivity for DEPAC relevant gases
3572	!-- Use default value
3573	diffusivity = 0.11e-4
3574	!
3575	!-- overwrite with known coefficients of diffusivity from Massman (1998)
3576	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
3577	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
3578	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
3579	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
3580	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3581	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3582	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3583	!
3584	!-- Get quasi-laminar boundary layer resistance rb:
3585	CALL get_rb_cell( (luv_dep == ilu_water_sea) .OR. (luv_dep == ilu_water_inland), &
3586	z0h_surf, ustar_surf, diffusivity, &
3587	rb )
3588	!
3589	!-- Get rc_tot
3590	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
3591	rh_surf, lai, sai, nwet, luv_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3592	r_aero_surf , rb )
3593	!
3594	!-- Calculate budget
3595	IF ( rc_tot <= 0.0 ) THEN
3596
3597	bud_luv(lsp) = 0.0_wp
3598
3599	ELSE
3600
3601	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3602
3603	bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3604	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3605	ENDIF
3606
3607	ENDIF
3608	ENDDO
3609	ENDIF
3610	!
3611	!-- Pavement
3612	IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 ) THEN
3613	!
3614	!-- No vegetation on pavements:
3615	lai = 0.0_wp
3616	sai = 0.0_wp
3617
3618	slinnfac = 1.0_wp
3619	!
3620	!-- Get vd
3621	DO lsp = 1, nvar
3622	!
3623	!-- Initialize
3624	vs = 0.0_wp
3625	vd_lu = 0.0_wp
3626	rs = 0.0_wp
3627	rb = 0.0_wp
3628	rc_tot = 0.0_wp
3629	IF ( spc_names(lsp) == 'PM10' ) THEN
3630	part_type = 1
3631	!
3632	!-- Sedimentation velocity:
3633	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3634	particle_pars(ind_p_size, part_type), &
3635	particle_pars(ind_p_slip, part_type), &
3636	visc)
3637
3638	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3639	vs, &
3640	particle_pars(ind_p_size, part_type), &
3641	particle_pars(ind_p_slip, part_type), &
3642	nwet, temp_tmp, dens, visc, &
3643	lup_dep, &
3644	r_aero_surf, ustar_surf )
3645
3646	bud_lup(lsp) = - conc_ijk(lsp) * &
3647	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3648
3649
3650	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
3651	part_type = 2
3652	!
3653	!-- Sedimentation velocity:
3654	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3655	particle_pars(ind_p_size, part_type), &
3656	particle_pars(ind_p_slip, part_type), &
3657	visc)
3658
3659	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3660	vs, &
3661	particle_pars(ind_p_size, part_type), &
3662	particle_pars(ind_p_slip, part_type), &
3663	nwet, temp_tmp, dens, visc, &
3664	lup_dep, &
3665	r_aero_surf, ustar_surf )
3666
3667	bud_lup(lsp) = - conc_ijk(lsp) * &
3668	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3669
3670	ELSE !<GASES
3671	!
3672	!-- Read spc_name of current species for gas parameter
3673
3674	IF ( ANY(pspecnames(:) == spc_names(lsp) ) ) THEN
3675	i_pspec = 0
3676	DO pspec = 1, nposp
3677	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
3678	i_pspec = pspec
3679	END IF
3680	ENDDO
3681
3682	ELSE
3683	!
3684	!-- For now species not deposited
3685	CYCLE
3686	ENDIF
3687	!
3688	!-- Factor used for conversion from ppb to ug/m3 :
3689	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3690	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3691	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
3692	!-- thus:
3693	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3694	!-- Use density at lowest layer:
3695
3696	ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3
3697	!
3698	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
3699	! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
3700	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
3701	!
3702	!-- Diffusivity for DEPAC relevant gases
3703	!-- Use default value
3704	diffusivity = 0.11e-4
3705	!
3706	!-- overwrite with known coefficients of diffusivity from Massman (1998)
3707	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
3708	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
3709	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
3710	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
3711	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3712	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3713	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3714	!
3715	!-- Get quasi-laminar boundary layer resistance rb:
3716	CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), &
3717	z0h_surf, ustar_surf, diffusivity, rb )
3718	!
3719	!-- Get rc_tot
3720	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, &
3721	solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2, &
3722	rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3723	r_aero_surf , rb )
3724	!
3725	!-- Calculate budget
3726	IF ( rc_tot <= 0.0 ) THEN
3727	bud_lup(lsp) = 0.0_wp
3728	ELSE
3729	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3730	bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3731	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3732	ENDIF
3733
3734	ENDIF
3735	ENDDO
3736	ENDIF
3737	!
3738	!-- Water
3739	IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 ) THEN
3740	!
3741	!-- No vegetation on water:
3742	lai = 0.0_wp
3743	sai = 0.0_wp
3744	slinnfac = 1.0_wp
3745	!
3746	!-- Get vd
3747	DO lsp = 1, nvar
3748	!
3749	!-- Initialize
3750	vs = 0.0_wp
3751	vd_lu = 0.0_wp
3752	rs = 0.0_wp
3753	rb = 0.0_wp
3754	rc_tot = 0.0_wp
3755	IF ( spc_names(lsp) == 'PM10' ) THEN
3756	part_type = 1
3757	!
3758	!-- Sedimentation velocity:
3759	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3760	particle_pars(ind_p_size, part_type), &
3761	particle_pars(ind_p_slip, part_type), &
3762	visc)
3763
3764	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3765	vs, &
3766	particle_pars(ind_p_size, part_type), &
3767	particle_pars(ind_p_slip, part_type), &
3768	nwet, temp_tmp, dens, visc, &
3769	luw_dep, &
3770	r_aero_surf, ustar_surf )
3771
3772	bud_luw(lsp) = - conc_ijk(lsp) * &
3773	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3774
3775	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
3776	part_type = 2
3777	!
3778	!-- Sedimentation velocity:
3779	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3780	particle_pars(ind_p_size, part_type), &
3781	particle_pars(ind_p_slip, part_type), &
3782	visc)
3783
3784	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3785	vs, &
3786	particle_pars(ind_p_size, part_type), &
3787	particle_pars(ind_p_slip, part_type), &
3788	nwet, temp_tmp, dens, visc, &
3789	luw_dep, &
3790	r_aero_surf, ustar_surf )
3791
3792	bud_luw(lsp) = - conc_ijk(lsp) * &
3793	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3794
3795	ELSE !<GASES
3796	!
3797	!-- Read spc_name of current species for gas PARAMETER
3798
3799	IF ( ANY(pspecnames(:) == spc_names(lsp) ) ) THEN
3800	i_pspec = 0
3801	DO pspec = 1, nposp
3802	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
3803	i_pspec = pspec
3804	END IF
3805	ENDDO
3806
3807	ELSE
3808	!
3809	!-- For now species not deposited
3810	CYCLE
3811	ENDIF
3812	!
3813	!-- Factor used for conversion from ppb to ug/m3 :
3814	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3815	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3816	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
3817	!-- thus:
3818	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3819	!-- Use density at lowest layer:
3820
3821	ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3
3822	!
3823	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
3824	!-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
3825	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
3826	!
3827	!-- Diffusivity for DEPAC relevant gases
3828	!-- Use default value
3829	diffusivity = 0.11e-4
3830	!
3831	!-- overwrite with known coefficients of diffusivity from Massman (1998)
3832	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
3833	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
3834	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
3835	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
3836	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3837	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3838	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3839	!
3840	!-- Get quasi-laminar boundary layer resistance rb:
3841	CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), &
3842	z0h_surf, ustar_surf, diffusivity, rb )
3843
3844	!-- Get rc_tot
3845	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, &
3846	solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2, &
3847	rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3848	r_aero_surf , rb )
3849	!
3850	!-- Calculate budget
3851	IF ( rc_tot <= 0.0 ) THEN
3852
3853	bud_luw(lsp) = 0.0_wp
3854
3855	ELSE
3856
3857	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3858
3859	bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3860	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
3861	ENDIF
3862
3863	ENDIF
3864	ENDDO
3865	ENDIF
3866
3867
3868	bud = 0.0_wp
3869	!
3870	!-- Calculate overall budget for surface m and adapt concentration
3871	DO lsp = 1, nspec
3872
3873	bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
3874	surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
3875	surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
3876	!
3877	!-- Compute new concentration:
3878	conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
3879
3880	chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
3881
3882	ENDDO
3883
3884	ENDIF
3885	!
3886	!-- For USM surfaces
3887
3888	IF ( match_usm ) THEN
3889	!
3890	!-- Get surface element information at i,j:
3891	m = surf_usm_h%start_index(j,i)
3892	k = surf_usm_h%k(m)
3893	!
3894	!-- Get needed variables for surface element m
3895	ustar_surf = surf_usm_h%us(m)
3896	z0h_surf = surf_usm_h%z0h(m)
3897	r_aero_surf = surf_usm_h%r_a(m)
3898	solar_rad = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
3899	lai = surf_usm_h%lai(m)
3900	sai = lai + 1
3901	!
3902	!-- For small grid spacing neglect R_a
3903	IF ( dzw(k) <= 1.0 ) THEN
3904	r_aero_surf = 0.0_wp
3905	ENDIF
3906	!
3907	!-- Initialize lu's
3908	luu_palm = 0
3909	luu_dep = 0
3910	lug_palm = 0
3911	lug_dep = 0
3912	lud_palm = 0
3913	lud_dep = 0
3914	!
3915	!-- Initialize budgets
3916	bud_luu = 0.0_wp
3917	bud_lug = 0.0_wp
3918	bud_lud = 0.0_wp
3919	!
3920	!-- Get land use for i,j and assign to DEPAC lu
3921	IF ( surf_usm_h%frac(ind_pav_green,m) > 0 ) THEN
3922	!
3923	!-- For green urban surfaces (e.g. green roofs
3924	!-- assume LU short grass
3925	lug_palm = ind_luv_s_grass
3926	IF ( lug_palm == ind_luv_user ) THEN
3927	message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3928	CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
3929	ELSEIF ( lug_palm == ind_luv_b_soil ) THEN
3930	lug_dep = 9
3931	ELSEIF ( lug_palm == ind_luv_mixed_crops ) THEN
3932	lug_dep = 2
3933	ELSEIF ( lug_palm == ind_luv_s_grass ) THEN
3934	lug_dep = 1
3935	ELSEIF ( lug_palm == ind_luv_ev_needle_trees ) THEN
3936	lug_dep = 4
3937	ELSEIF ( lug_palm == ind_luv_de_needle_trees ) THEN
3938	lug_dep = 4
3939	ELSEIF ( lug_palm == ind_luv_ev_broad_trees ) THEN
3940	lug_dep = 12
3941	ELSEIF ( lug_palm == ind_luv_de_broad_trees ) THEN
3942	lug_dep = 5
3943	ELSEIF ( lug_palm == ind_luv_t_grass ) THEN
3944	lug_dep = 1
3945	ELSEIF ( lug_palm == ind_luv_desert ) THEN
3946	lug_dep = 9
3947	ELSEIF ( lug_palm == ind_luv_tundra ) THEN
3948	lug_dep = 8
3949	ELSEIF ( lug_palm == ind_luv_irr_crops ) THEN
3950	lug_dep = 2
3951	ELSEIF ( lug_palm == ind_luv_semidesert ) THEN
3952	lug_dep = 8
3953	ELSEIF ( lug_palm == ind_luv_ice ) THEN
3954	lug_dep = 10
3955	ELSEIF ( lug_palm == ind_luv_marsh ) THEN
3956	lug_dep = 8
3957	ELSEIF ( lug_palm == ind_luv_ev_shrubs ) THEN
3958	lug_dep = 14
3959	ELSEIF ( lug_palm == ind_luv_de_shrubs ) THEN
3960	lug_dep = 14
3961	ELSEIF ( lug_palm == ind_luv_mixed_forest ) THEN
3962	lug_dep = 4
3963	ELSEIF ( lug_palm == ind_luv_intrup_forest ) THEN
3964	lug_dep = 8
3965	ENDIF
3966	ENDIF
3967
3968	IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 ) THEN
3969	!
3970	!-- For walls in USM assume concrete walls/roofs,
3971	!-- assumed LU class desert as also assumed for
3972	!-- pavements in LSM
3973	luu_palm = ind_lup_conc
3974	IF ( luu_palm == ind_lup_user ) THEN
3975	message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3976	CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
3977	ELSEIF ( luu_palm == ind_lup_asph_conc ) THEN
3978	luu_dep = 9
3979	ELSEIF ( luu_palm == ind_lup_asph ) THEN
3980	luu_dep = 9
3981	ELSEIF ( luu_palm == ind_lup_conc ) THEN
3982	luu_dep = 9
3983	ELSEIF ( luu_palm == ind_lup_sett ) THEN
3984	luu_dep = 9
3985	ELSEIF ( luu_palm == ind_lup_pav_stones ) THEN
3986	luu_dep = 9
3987	ELSEIF ( luu_palm == ind_lup_cobblest ) THEN
3988	luu_dep = 9
3989	ELSEIF ( luu_palm == ind_lup_metal ) THEN
3990	luu_dep = 9
3991	ELSEIF ( luu_palm == ind_lup_wood ) THEN
3992	luu_dep = 9
3993	ELSEIF ( luu_palm == ind_lup_gravel ) THEN
3994	luu_dep = 9
3995	ELSEIF ( luu_palm == ind_lup_f_gravel ) THEN
3996	luu_dep = 9
3997	ELSEIF ( luu_palm == ind_lup_pebblest ) THEN
3998	luu_dep = 9
3999	ELSEIF ( luu_palm == ind_lup_woodchips ) THEN
4000	luu_dep = 9
4001	ELSEIF ( luu_palm == ind_lup_tartan ) THEN
4002	luu_dep = 9
4003	ELSEIF ( luu_palm == ind_lup_art_turf ) THEN
4004	luu_dep = 9
4005	ELSEIF ( luu_palm == ind_lup_clay ) THEN
4006	luu_dep = 9
4007	ENDIF
4008	ENDIF
4009
4010	IF ( surf_usm_h%frac(ind_wat_win,m) > 0 ) THEN
4011	!
4012	!-- For windows in USM assume metal as this is
4013	!-- as close as we get, assumed LU class desert
4014	!-- as also assumed for pavements in LSM
4015	lud_palm = ind_lup_metal
4016	IF ( lud_palm == ind_lup_user ) THEN
4017	message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
4018	CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
4019	ELSEIF ( lud_palm == ind_lup_asph_conc ) THEN
4020	lud_dep = 9
4021	ELSEIF ( lud_palm == ind_lup_asph ) THEN
4022	lud_dep = 9
4023	ELSEIF ( lud_palm == ind_lup_conc ) THEN
4024	lud_dep = 9
4025	ELSEIF ( lud_palm == ind_lup_sett ) THEN
4026	lud_dep = 9
4027	ELSEIF ( lud_palm == ind_lup_pav_stones ) THEN
4028	lud_dep = 9
4029	ELSEIF ( lud_palm == ind_lup_cobblest ) THEN
4030	lud_dep = 9
4031	ELSEIF ( lud_palm == ind_lup_metal ) THEN
4032	lud_dep = 9
4033	ELSEIF ( lud_palm == ind_lup_wood ) THEN
4034	lud_dep = 9
4035	ELSEIF ( lud_palm == ind_lup_gravel ) THEN
4036	lud_dep = 9
4037	ELSEIF ( lud_palm == ind_lup_f_gravel ) THEN
4038	lud_dep = 9
4039	ELSEIF ( lud_palm == ind_lup_pebblest ) THEN
4040	lud_dep = 9
4041	ELSEIF ( lud_palm == ind_lup_woodchips ) THEN
4042	lud_dep = 9
4043	ELSEIF ( lud_palm == ind_lup_tartan ) THEN
4044	lud_dep = 9
4045	ELSEIF ( lud_palm == ind_lup_art_turf ) THEN
4046	lud_dep = 9
4047	ELSEIF ( lud_palm == ind_lup_clay ) THEN
4048	lud_dep = 9
4049	ENDIF
4050	ENDIF
4051	!
4052	!-- @TODO: Activate these lines as soon as new ebsolver branch is merged:
4053	!-- Set wetness indicator to dry or wet for usm vegetation or pavement
4054	!IF ( surf_usm_h%c_liq(m) > 0 ) THEN
4055	! nwet = 1
4056	!ELSE
4057	nwet = 0
4058	!ENDIF
4059	!
4060	!-- Compute length of time step
4061	IF ( call_chem_at_all_substeps ) THEN
4062	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
4063	ELSE
4064	dt_chem = dt_3d
4065	ENDIF
4066
4067	dh = dzw(k)
4068	inv_dh = 1.0_wp / dh
4069	dt_dh = dt_chem/dh
4070	!
4071	!-- Concentration at i,j,k
4072	DO lsp = 1, nspec
4073	conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
4074	ENDDO
4075	!
4076	!-- Temperature at i,j,k
4077	temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
4078
4079	ts = temp_tmp - 273.15 !< in degrees celcius
4080	!
4081	!-- Viscosity of air
4082	visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
4083	!
4084	!-- Air density at k
4085	dens = rho_air_zw(k)
4086	!
4087	!-- Calculate relative humidity from specific humidity for DEPAC
4088	qv_tmp = MAX(q(k,j,i),0.0_wp)
4089	rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
4090	!
4091	!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
4092	!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
4093	!
4094	!-- Walls/roofs
4095	IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 ) THEN
4096	!
4097	!-- No vegetation on non-green walls:
4098	lai = 0.0_wp
4099	sai = 0.0_wp
4100
4101	slinnfac = 1.0_wp
4102	!
4103	!-- Get vd
4104	DO lsp = 1, nvar
4105	!
4106	!-- Initialize
4107	vs = 0.0_wp
4108	vd_lu = 0.0_wp
4109	rs = 0.0_wp
4110	rb = 0.0_wp
4111	rc_tot = 0.0_wp
4112	IF (spc_names(lsp) == 'PM10' ) THEN
4113	part_type = 1
4114	!
4115	!-- Sedimentation velocity
4116	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4117	particle_pars(ind_p_size, part_type), &
4118	particle_pars(ind_p_slip, part_type), &
4119	visc)
4120
4121	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4122	vs, &
4123	particle_pars(ind_p_size, part_type), &
4124	particle_pars(ind_p_slip, part_type), &
4125	nwet, temp_tmp, dens, visc, &
4126	luu_dep, &
4127	r_aero_surf, ustar_surf )
4128
4129	bud_luu(lsp) = - conc_ijk(lsp) * &
4130	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4131
4132	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
4133	part_type = 2
4134	!
4135	!-- Sedimentation velocity
4136	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4137	particle_pars(ind_p_size, part_type), &
4138	particle_pars(ind_p_slip, part_type), &
4139	visc)
4140
4141	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4142	vs, &
4143	particle_pars(ind_p_size, part_type), &
4144	particle_pars(ind_p_slip, part_type), &
4145	nwet, temp_tmp, dens, visc, &
4146	luu_dep , &
4147	r_aero_surf, ustar_surf )
4148
4149	bud_luu(lsp) = - conc_ijk(lsp) * &
4150	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4151
4152	ELSE !< GASES
4153	!
4154	!-- Read spc_name of current species for gas parameter
4155
4156	IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN
4157	i_pspec = 0
4158	DO pspec = 1, nposp
4159	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
4160	i_pspec = pspec
4161	END IF
4162	ENDDO
4163	ELSE
4164	!
4165	!-- For now species not deposited
4166	CYCLE
4167	ENDIF
4168	!
4169	!-- Factor used for conversion from ppb to ug/m3 :
4170	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4171	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4172	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
4173	!-- thus:
4174	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4175	!-- Use density at k:
4176
4177	ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3
4178
4179	!
4180	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
4181	!-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
4182	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
4183	!
4184	!-- Diffusivity for DEPAC relevant gases
4185	!-- Use default value
4186	diffusivity = 0.11e-4
4187	!
4188	!-- overwrite with known coefficients of diffusivity from Massman (1998)
4189	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4190	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
4191	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
4192	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
4193	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4194	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4195	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4196	!
4197	!-- Get quasi-laminar boundary layer resistance rb:
4198	CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), &
4199	z0h_surf, ustar_surf, diffusivity, &
4200	rb )
4201	!
4202	!-- Get rc_tot
4203	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, &
4204	solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2, &
4205	rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
4206	r_aero_surf, rb )
4207	!
4208	!-- Calculate budget
4209	IF ( rc_tot <= 0.0 ) THEN
4210
4211	bud_luu(lsp) = 0.0_wp
4212
4213	ELSE
4214
4215	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4216
4217	bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4218	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4219	ENDIF
4220
4221	ENDIF
4222	ENDDO
4223	ENDIF
4224	!
4225	!-- Green usm surfaces
4226	IF ( surf_usm_h%frac(ind_pav_green,m) > 0 ) THEN
4227
4228	!
4229	!-- No vegetation on bare soil, desert or ice:
4230	IF ( ( lug_palm == ind_luv_b_soil ) .OR. &
4231	( lug_palm == ind_luv_desert ) .OR. &
4232	( lug_palm == ind_luv_ice ) ) THEN
4233
4234	lai = 0.0_wp
4235	sai = 0.0_wp
4236
4237	ENDIF
4238
4239
4240	slinnfac = 1.0_wp
4241	!
4242	!-- Get vd
4243	DO lsp = 1, nvar
4244	!
4245	!-- Initialize
4246	vs = 0.0_wp
4247	vd_lu = 0.0_wp
4248	rs = 0.0_wp
4249	rb = 0.0_wp
4250	rc_tot = 0.0_wp
4251	IF ( spc_names(lsp) == 'PM10' ) THEN
4252	part_type = 1
4253	!
4254	!-- Sedimentation velocity
4255	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4256	particle_pars(ind_p_size, part_type), &
4257	particle_pars(ind_p_slip, part_type), &
4258	visc)
4259
4260	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4261	vs, &
4262	particle_pars(ind_p_size, part_type), &
4263	particle_pars(ind_p_slip, part_type), &
4264	nwet, temp_tmp, dens, visc, &
4265	lug_dep, &
4266	r_aero_surf, ustar_surf )
4267
4268	bud_lug(lsp) = - conc_ijk(lsp) * &
4269	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4270
4271	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
4272	part_type = 2
4273	!
4274	!-- Sedimentation velocity
4275	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4276	particle_pars(ind_p_size, part_type), &
4277	particle_pars(ind_p_slip, part_type), &
4278	visc)
4279
4280	CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4281	vs, &
4282	particle_pars(ind_p_size, part_type), &
4283	particle_pars(ind_p_slip, part_type), &
4284	nwet, temp_tmp, dens, visc, &
4285	lug_dep, &
4286	r_aero_surf, ustar_surf )
4287
4288	bud_lug(lsp) = - conc_ijk(lsp) * &
4289	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4290
4291	ELSE !< GASES
4292	!
4293	!-- Read spc_name of current species for gas parameter
4294
4295	IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN
4296	i_pspec = 0
4297	DO pspec = 1, nposp
4298	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
4299	i_pspec = pspec
4300	END IF
4301	ENDDO
4302	ELSE
4303	!
4304	!-- For now species not deposited
4305	CYCLE
4306	ENDIF
4307	!
4308	!-- Factor used for conversion from ppb to ug/m3 :
4309	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4310	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4311	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
4312	!-- thus:
4313	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4314	!-- Use density at k:
4315
4316	ppm2ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3
4317	!
4318	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
4319	! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
4320	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
4321	!
4322	!-- Diffusivity for DEPAC relevant gases
4323	!-- Use default value
4324	diffusivity = 0.11e-4
4325	!
4326	!-- overwrite with known coefficients of diffusivity from Massman (1998)
4327	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4328	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
4329	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
4330	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
4331	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4332	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4333	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4334	!
4335	!-- Get quasi-laminar boundary layer resistance rb:
4336	CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), &
4337	z0h_surf, ustar_surf, diffusivity, &
4338	rb )
4339	!
4340	!-- Get rc_tot
4341	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, &
4342	solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2, &
4343	rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
4344	r_aero_surf , rb )
4345	!
4346	!-- Calculate budget
4347	IF ( rc_tot <= 0.0 ) THEN
4348
4349	bud_lug(lsp) = 0.0_wp
4350
4351	ELSE
4352
4353	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4354
4355	bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4356	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4357	ENDIF
4358
4359	ENDIF
4360	ENDDO
4361	ENDIF
4362	!
4363	!-- Windows
4364	IF ( surf_usm_h%frac(ind_wat_win,m) > 0 ) THEN
4365	!
4366	!-- No vegetation on windows:
4367	lai = 0.0_wp
4368	sai = 0.0_wp
4369
4370	slinnfac = 1.0_wp
4371	!
4372	!-- Get vd
4373	DO lsp = 1, nvar
4374	!
4375	!-- Initialize
4376	vs = 0.0_wp
4377	vd_lu = 0.0_wp
4378	rs = 0.0_wp
4379	rb = 0.0_wp
4380	rc_tot = 0.0_wp
4381	IF ( spc_names(lsp) == 'PM10' ) THEN
4382	part_type = 1
4383	!
4384	!-- Sedimentation velocity
4385	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4386	particle_pars(ind_p_size, part_type), &
4387	particle_pars(ind_p_slip, part_type), &
4388	visc)
4389
4390	CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4391	particle_pars(ind_p_size, part_type), &
4392	particle_pars(ind_p_slip, part_type), &
4393	nwet, temp_tmp, dens, visc, &
4394	lud_dep, r_aero_surf, ustar_surf )
4395
4396	bud_lud(lsp) = - conc_ijk(lsp) * &
4397	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4398
4399	ELSEIF ( spc_names(lsp) == 'PM25' ) THEN
4400	part_type = 2
4401	!
4402	!-- Sedimentation velocity
4403	vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4404	particle_pars(ind_p_size, part_type), &
4405	particle_pars(ind_p_slip, part_type), &
4406	visc)
4407
4408	CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4409	particle_pars(ind_p_size, part_type), &
4410	particle_pars(ind_p_slip, part_type), &
4411	nwet, temp_tmp, dens, visc, &
4412	lud_dep, &
4413	r_aero_surf, ustar_surf )
4414
4415	bud_lud(lsp) = - conc_ijk(lsp) * &
4416	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4417
4418	ELSE !< GASES
4419	!
4420	!-- Read spc_name of current species for gas PARAMETER
4421
4422	IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN
4423	i_pspec = 0
4424	DO pspec = 1, nposp
4425	IF ( pspecnames(pspec) == spc_names(lsp) ) THEN
4426	i_pspec = pspec
4427	END IF
4428	ENDDO
4429	ELSE
4430	!
4431	!-- For now species not deposited
4432	CYCLE
4433	ENDIF
4434	!
4435	!-- Factor used for conversion from ppb to ug/m3 :
4436	!-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4437	!-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4438	!-- c 1e-9 xm_tracer 1e9 / xm_air dens
4439	!-- thus:
4440	!-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4441	!-- Use density at k:
4442
4443	ppm2ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3
4444	!
4445	!-- Atmospheric concentration in DEPAC is requested in ug/m3:
4446	!-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole
4447	conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3
4448	!
4449	!-- Diffusivity for DEPAC relevant gases
4450	!-- Use default value
4451	diffusivity = 0.11e-4
4452	!
4453	!-- overwrite with known coefficients of diffusivity from Massman (1998)
4454	IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4455	IF ( spc_names(lsp) == 'NO' ) diffusivity = 0.199e-4
4456	IF ( spc_names(lsp) == 'O3' ) diffusivity = 0.144e-4
4457	IF ( spc_names(lsp) == 'CO' ) diffusivity = 0.176e-4
4458	IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4459	IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4460	IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4461	!
4462	!-- Get quasi-laminar boundary layer resistance rb:
4463	CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), &
4464	z0h_surf, ustar_surf, diffusivity, rb )
4465	!
4466	!-- Get rc_tot
4467	CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, &
4468	solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2, &
4469	rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
4470	r_aero_surf , rb )
4471	!
4472	!-- Calculate budget
4473	IF ( rc_tot <= 0.0 ) THEN
4474
4475	bud_lud(lsp) = 0.0_wp
4476
4477	ELSE
4478
4479	vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4480
4481	bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4482	(1.0_wp - exp(-vd_lu * dt_dh )) * dh
4483	ENDIF
4484
4485	ENDIF
4486	ENDDO
4487	ENDIF
4488
4489
4490	bud = 0.0_wp
4491	!
4492	!-- Calculate overall budget for surface m and adapt concentration
4493	DO lsp = 1, nspec
4494
4495
4496	bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
4497	surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
4498	surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
4499	!
4500	!-- Compute new concentration
4501	conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4502
4503	chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
4504
4505	ENDDO
4506
4507	ENDIF
4508
4509
4510	END SUBROUTINE chem_depo
4511
4512
4513	!------------------------------------------------------------------------------!
4514	! Description:
4515	! ------------
4516	!> Subroutine to compute total canopy (or surface) resistance Rc for gases
4517	!>
4518	!> DEPAC:
4519	!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
4520	!> module of the LOTOS-EUROS model (Manders et al., 2017)
4521	!>
4522	!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
4523	!> van Zanten et al., 2010.
4524	!------------------------------------------------------------------------------!
4525	SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi, &
4526	rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity, &
4527	ra, rb )
4528	!
4529	!-- Some of depac arguments are OPTIONAL:
4530	!-- A. compute Rc_tot without compensation points (ccomp_tot will be zero):
4531	!-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
4532	!-- B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
4533	!-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4534	!-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
4535	!--
4536	!-- C. compute effective Rc based on compensation points (used for OPS):
4537	!-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4538	!-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4539	!-- ra, rb, rc_eff)
4540	!-- X1. Extra (OPTIONAL) output variables:
4541	!-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4542	!-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4543	!-- ra, rb, rc_eff, &
4544	!-- gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
4545	!-- X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
4546	!-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4547	!-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4548	!-- ra, rb, rc_eff, &
4549	!-- gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
4550	!-- calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
4551
4552
4553	CHARACTER(LEN=*), INTENT(IN) :: compnam !< component name
4554	!< 'HNO3','NO','NO2','O3','SO2','NH3'
4555	INTEGER(iwp), INTENT(IN) :: day_of_year !< day of year, 1 ... 365 (366)
4556	INTEGER(iwp), INTENT(IN) :: nwet !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4557	INTEGER(iwp), INTENT(IN) :: lu !< land use type, lu = 1,...,nlu
4558	INTEGER(iwp), INTENT(IN) :: iratns !< index for NH3/SO2 ratio used for SO2:
4559	!< iratns = 1: low NH3/SO2
4560	!< iratns = 2: high NH3/SO2
4561	!< iratns = 3: very low NH3/SO2
4562	REAL(wp), INTENT(IN) :: lat !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
4563	REAL(wp), INTENT(IN) :: t !< temperature (C)
4564	REAL(wp), INTENT(IN) :: ust !< friction velocity (m/s)
4565	REAL(wp), INTENT(IN) :: solar_rad !< solar radiation, dirict+diffuse (W/m2)
4566	REAL(wp), INTENT(IN) :: sinphi !< sin of solar elevation angle
4567	REAL(wp), INTENT(IN) :: rh !< relative humidity (%)
4568	REAL(wp), INTENT(IN) :: lai !< one-sidedleaf area index (-)
4569	REAL(wp), INTENT(IN) :: sai !< surface area index (-) (lai + branches and stems)
4570	REAL(wp), INTENT(IN) :: diffusivity !< diffusivity
4571	REAL(wp), INTENT(IN) :: p !< pressure (Pa)
4572	REAL(wp), INTENT(IN) :: conc_ijk_ugm3 !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
4573	REAL(wp), INTENT(IN) :: ra !< aerodynamic resistance (s/m)
4574	REAL(wp), INTENT(IN) :: rb !< boundary layer resistance (s/m)
4575
4576	REAL(wp), INTENT(OUT) :: rc_tot !< total canopy resistance Rc (s/m)
4577	REAL(wp), INTENT(OUT) :: ccomp_tot !< total compensation point (ug/m3)
4578	! !< [= 0 for species that don't have a compensation
4579	!-- Local variables:
4580	!
4581	!-- Component number taken from component name, paramteres matched with include files
4582	INTEGER(iwp) :: icmp
4583	!
4584	!-- Component numbers:
4585	INTEGER(iwp), PARAMETER :: icmp_o3 = 1
4586	INTEGER(iwp), PARAMETER :: icmp_so2 = 2
4587	INTEGER(iwp), PARAMETER :: icmp_no2 = 3
4588	INTEGER(iwp), PARAMETER :: icmp_no = 4
4589	INTEGER(iwp), PARAMETER :: icmp_nh3 = 5
4590	INTEGER(iwp), PARAMETER :: icmp_co = 6
4591	INTEGER(iwp), PARAMETER :: icmp_no3 = 7
4592	INTEGER(iwp), PARAMETER :: icmp_hno3 = 8
4593	INTEGER(iwp), PARAMETER :: icmp_n2o5 = 9
4594	INTEGER(iwp), PARAMETER :: icmp_h2o2 = 10
4595
4596	LOGICAL :: ready !< Rc has been set:
4597	!< = 1 -> constant Rc
4598	!< = 2 -> temperature dependent Rc
4599	!
4600	!-- Vegetation indicators:
4601	LOGICAL :: lai_present !< leaves are present for current land use type
4602