Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |