Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |