Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |