Home

Context Navigation

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

Last change on this file since 3638 was 3638, checked in by forkel, 6 years ago
chemistry_model_mod: Added missing conversion factor fr2ppm for qvap
Property svn:keywords set to `Id`
File size: 216.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |