Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |