Home

Context Navigation

chemistry_model_mod.f90 @ 4004

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