Home

Context Navigation

source: palm/trunk/SOURCE/chemistry_model_mod.f90 @ 3887

Last change on this file since 3887 was 3887, checked in by schwenkel, 6 years ago
bugfix for chemistry_model_mod via introducing module_interface_exchange_horiz
Property svn:keywords set to `Id`
File size: 225.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |