Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |