Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |