Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |