Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |