Home

Context Navigation

chemistry_model_mod.f90 @ 4167

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