Home

Context Navigation

chemistry_model_mod.f90 @ 4029

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