Home

Context Navigation

chemistry_model_mod.f90 @ 4102

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