Home

Context Navigation

chemistry_model_mod.f90 @ 4306

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