Home

Context Navigation

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

Last change on this file since 3543 was 3543, checked in by suehring, 6 years ago
changes from last commit documented
Property svn:keywords set to `Id`
File size: 224.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |