Home

Context Navigation

chemistry_model_mod.f90 @ 3543

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