Home

Context Navigation

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

Last change on this file since 3885 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 224.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |