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

Last change on this file since 4110 was 4110, checked in by suehring, 4 years ago

last changes documented

  • Property svn:keywords set to Id
File size: 239.1 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-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Karlsruhe Institute of Technology
19! Copyright 2017-2019 Freie Universitaet Berlin
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! -----------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chemistry_model_mod.f90 4110 2019-07-22 17:05:21Z suehring $
29! - Decycling boundary conditions are only set at the ghost points not on the
30!   prognostic grid points
31! - Allocation and initialization of special advection flags cs_advc_flags_s
32!   used for chemistry. These are exclusively used for chemical species in
33!   order to distinguish from the usually-used flags which might be different
34!   when decycling is applied in combination with cyclic boundary conditions.
35!   Moreover, cs_advc_flags_s considers extended zones around buildings where
36!   first-order upwind scheme is applied for the horizontal advection terms,
37!   in order to overcome high concentration peaks due to stationary numerical
38!   oscillations caused by horizontal advection discretization.
39!
40! 4109 2019-07-22 17:00:34Z suehring
41! Slightly revise setting of boundary conditions at horizontal walls, use
42! data-structure offset index instead of pre-calculate it for each facing
43!
44! 4080 2019-07-09 18:17:37Z suehring
45! Restore accidantly removed limitation to positive values
46!
47! 4079 2019-07-09 18:04:41Z suehring
48! Application of monotonic flux limiter for the vertical scalar advection
49! up to the topography top (only for the cache-optimized version at the
50! moment).
51!
52! 4069 2019-07-01 14:05:51Z Giersch
53! Masked output running index mid has been introduced as a local variable to
54! avoid runtime error (Loop variable has been modified) in time_integration
55!
56! 4029 2019-06-14 14:04:35Z raasch
57! nest_chemistry option removed
58!
59! 4004 2019-05-24 11:32:38Z suehring
60! in subroutine chem_parin check emiss_lod / mod_emis only
61! when emissions_anthropogenic is activated in namelist (E.C. Chan)
62!
63! 3968 2019-05-13 11:04:01Z suehring
64! - added "emiss_lod" which serves the same function as "mode_emis"
65!   both will be synchronized with emiss_lod having pirority over
66!   mode_emis (see informational messages)
67! - modified existing error message and introduced new informational messages
68!    - CM0436 - now also applies to invalid LOD settings
69!    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
70!    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
71!
72! 3930 2019-04-24 14:57:18Z forkel
73! Changed chem_non_transport_physics to chem_non_advective_processes
74!
75!
76! 3929 2019-04-24 12:52:08Z banzhafs
77! Correct/complete module_interface introduction for chemistry model
78! Add subroutine chem_exchange_horiz_bounds
79! Bug fix deposition
80!
81! 3784 2019-03-05 14:16:20Z banzhafs
82! 2D output of emission fluxes
83!
84! 3784 2019-03-05 14:16:20Z banzhafs
85! Bugfix, uncomment erroneous commented variable used for dry deposition.
86! Bugfix in 3D emission output.
87!
88! 3784 2019-03-05 14:16:20Z banzhafs
89! Changes related to global restructuring of location messages and introduction
90! of additional debug messages
91!
92! 3784 2019-03-05 14:16:20Z banzhafs
93! some formatting of the deposition code
94!
95! 3784 2019-03-05 14:16:20Z banzhafs
96! some formatting
97!
98! 3784 2019-03-05 14:16:20Z banzhafs
99! added cs_mech to USE chem_gasphase_mod
100!
101! 3784 2019-03-05 14:16:20Z banzhafs
102! renamed get_mechanismname to get_mechanism_name
103! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
104!
105! 3784 2019-03-05 14:16:20Z banzhafs
106! Unused variables removed/taken care of
107!
108!
109! 3784 2019-03-05 14:16:20Z forkel
110! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
111!
112!
113! 3783 2019-03-05 13:23:50Z forkel
114! Removed forgotte write statements an some unused variables (did not touch the
115! parts related to deposition)
116!
117!
118! 3780 2019-03-05 11:19:45Z forkel
119! Removed READ from unit 10, added CALL get_mechanismname
120!
121!
122! 3767 2019-02-27 08:18:02Z raasch
123! unused variable for file index removed from rrd-subroutines parameter list
124!
125! 3738 2019-02-12 17:00:45Z suehring
126! Clean-up debug prints
127!
128! 3737 2019-02-12 16:57:06Z suehring
129! Enable mesoscale offline nesting for chemistry variables as well as
130! initialization of chemistry via dynamic input file.
131!
132! 3719 2019-02-06 13:10:18Z kanani
133! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call
134! to prognostic_equations for better overview
135!
136! 3700 2019-01-26 17:03:42Z knoop
137! Some interface calls moved to module_interface + cleanup
138!
139! 3664 2019-01-09 14:00:37Z forkel
140! Replaced misplaced location message by @todo
141!
142!
143! 3654 2019-01-07 16:31:57Z suehring
144! Disable misplaced location message
145!
146! 3652 2019-01-07 15:29:59Z forkel
147! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
148!
149!
150! 3646 2018-12-28 17:58:49Z kanani
151! Bugfix: use time_since_reference_point instead of simulated_time (relevant
152! when using wall/soil spinup)
153!
154! 3643 2018-12-24 13:16:19Z knoop
155! Bugfix: set found logical correct in chem_data_output_2d
156!
157! 3638 2018-12-20 13:18:23Z forkel
158! Added missing conversion factor fr2ppm for qvap
159!
160!
161! 3637 2018-12-20 01:51:36Z knoop
162! Implementation of the PALM module interface
163!
164! 3636 2018-12-19 13:48:34Z raasch
165! nopointer option removed
166!
167! 3611 2018-12-07 14:14:11Z banzhafs
168! Minor formatting             
169!
170! 3600 2018-12-04 13:49:07Z banzhafs
171! Code update to comply PALM coding rules           
172! Bug fix in par_dir_diff subroutine                 
173! Small fixes (corrected 'conastant', added 'Unused')
174!
175! 3586 2018-11-30 13:20:29Z dom_dwd_user
176! Changed character lenth of name in species_def and photols_def to 15
177!
178! 3570 2018-11-27 17:44:21Z kanani
179! resler:
180! Break lines at 132 characters
181!
182! 3543 2018-11-20 17:06:15Z suehring
183! Remove tabs
184!
185! 3542 2018-11-20 17:04:13Z suehring
186! working precision added to make code Fortran 2008 conform
187!
188! 3458 2018-10-30 14:51:23Z kanani
189! from chemistry branch r3443, banzhafs, basit:
190! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo
191!
192! bug fix in chem_depo: allow different surface fractions for one
193! surface element and set lai to zero for non vegetated surfaces
194! bug fixed in chem_data_output_2d
195! bug fix in chem_depo subroutine
196! added code on deposition of gases and particles
197! removed cs_profile_name from chem_parin
198! bug fixed in output profiles and code cleaned
199!
200! 3449 2018-10-29 19:36:56Z suehring
201! additional output - merged from branch resler
202!
203! 3438 2018-10-28 19:31:42Z pavelkrc
204! Add terrain-following masked output
205!
206! 3373 2018-10-18 15:25:56Z kanani
207! Remove MPI_Abort, replace by message
208!
209! 3318 2018-10-08 11:43:01Z sward
210! Fixed faulty syntax of message string
211!
212! 3298 2018-10-02 12:21:11Z kanani
213! Add remarks (kanani)
214! Merge with trunk, replaced cloud_physics by bulk_cloud_model         28.09.2018 forkel
215! Subroutines header and chem_check_parameters added                   25.09.2018 basit
216! Removed chem_emission routine now declared in chem_emissions.f90     30.07.2018 ERUSSO
217! Introduced emissions namelist parameters                             30.07.2018 ERUSSO
218!
219! Timestep steering added in subroutine chem_integrate_ij and
220! output of chosen solver in chem_parin added                          30.07.2018 ketelsen
221!
222! chem_check_data_output_pr: added unit for PM compounds               20.07.2018 forkel
223! replaced : by nzb+1:nzt for pt,q,ql (found by kk)                    18.07.2018 forkel
224! debugged restart run for chem species               06.07.2018 basit
225! reorganized subroutines in alphabetical order.      27.06.2018 basit
226! subroutine chem_parin updated for profile output    27.06.2018 basit
227! Added humidity arrays to USE section and tmp_qvap in chem_integrate  26.6.2018 forkel
228! Merged chemistry with with trunk (nzb_do and nzt_do in 3d output)    26.6.2018 forkel
229!
230! reorganized subroutines in alphabetical order.      basit 22.06.2018
231! subroutine chem_parin updated for profile output    basit 22.06.2018
232! subroutine chem_statistics added                 
233! subroutine chem_check_data_output_pr add            21.06.2018 basit
234! subroutine chem_data_output_mask added              20.05.2018 basit
235! subroutine chem_data_output_2d added                20.05.2018 basit
236! subroutine chem_statistics added                    04.06.2018 basit
237! subroutine chem_emissions: Set cssws to zero before setting values 20.03.2018 forkel
238! subroutine chem_emissions: Introduced different conversion factors
239! for PM and gaseous compounds                                    15.03.2018 forkel
240! subroutine chem_emissions updated to take variable number of chem_spcs and
241! emission factors.                                               13.03.2018 basit
242! chem_boundary_conds_decycle improved.                           05.03.2018 basit
243! chem_boundary_conds_decycle subroutine added                    21.02.2018 basit
244! chem_init_profiles subroutines re-activated after correction    21.02.2018 basit
245!
246!
247! 3293 2018-09-28 12:45:20Z forkel
248! Modularization of all bulk cloud physics code components
249!
250! 3248 2018-09-14 09:42:06Z sward
251! Minor formating changes
252!
253! 3246 2018-09-13 15:14:50Z sward
254! Added error handling for input namelist via parin_fail_message
255!
256! 3241 2018-09-12 15:02:00Z raasch
257! +nest_chemistry
258!
259! 3209 2018-08-27 16:58:37Z suehring
260! Rename flags indicating outflow boundary conditions
261!
262! 3182 2018-07-27 13:36:03Z suehring
263! Revise output of surface quantities in case of overhanging structures
264!
265! 3045 2018-05-28 07:55:41Z Giersch
266! error messages revised
267!
268! 3014 2018-05-09 08:42:38Z maronga
269! Bugfix: nzb_do and nzt_do were not used for 3d data output
270!
271! 3004 2018-04-27 12:33:25Z Giersch
272! Comment concerning averaged data output added
273!
274! 2932 2018-03-26 09:39:22Z maronga
275! renamed chemistry_par to chemistry_parameters
276!
277! 2894 2018-03-15 09:17:58Z Giersch
278! Calculations of the index range of the subdomain on file which overlaps with
279! the current subdomain are already done in read_restart_data_mod,
280! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
281! renamed to chem_rrd_local, chem_write_var_list was renamed to
282! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
283! chem_skip_var_list has been removed, variable named found has been
284! introduced for checking if restart data was found, reading of restart strings
285! has been moved completely to read_restart_data_mod, chem_rrd_local is already
286! inside the overlap loop programmed in read_restart_data_mod, todo list has
287! bees extended, redundant characters in chem_wrd_local have been removed,
288! the marker *** end chemistry *** is not necessary anymore, strings and their
289! respective lengths are written out and read now in case of restart runs to
290! get rid of prescribed character lengths
291!
292! 2815 2018-02-19 11:29:57Z suehring
293! Bugfix in restart mechanism,
294! rename chem_tendency to chem_prognostic_equations,
295! implement vector-optimized version of chem_prognostic_equations,
296! some clean up (incl. todo list)
297!
298! 2773 2018-01-30 14:12:54Z suehring
299! Declare variables required for nesting as public
300!
301! 2772 2018-01-29 13:10:35Z suehring
302! Bugfix in string handling
303!
304! 2768 2018-01-24 15:38:29Z kanani
305! Shorten lines to maximum length of 132 characters
306!
307! 2766 2018-01-22 17:17:47Z kanani
308! Removed preprocessor directive __chem
309!
310! 2756 2018-01-16 18:11:14Z suehring
311! Fill values in 3D output introduced.
312!
313! 2718 2018-01-02 08:49:38Z maronga
314! Initial revision
315!
316!
317!
318!
319! Authors:
320! --------
321! @author Renate Forkel
322! @author Farah Kanani-Suehring
323! @author Klaus Ketelsen
324! @author Basit Khan
325! @author Sabine Banzhaf
326!
327!
328!------------------------------------------------------------------------------!
329! Description:
330! ------------
331!> Chemistry model for PALM-4U
332!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
333!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
334!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
335!>       allowed to use the chemistry model in a precursor run and additionally
336!>       not using it in a main run
337!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
338!> @todo Separate boundary conditions for each chem spcs to be implemented
339!> @todo Currently only total concentration are calculated. Resolved, parameterized
340!>       and chemistry fluxes although partially and some completely coded but
341!>       are not operational/activated in this version. bK.
342!> @todo slight differences in passive scalar and chem spcs when chem reactions
343!>       turned off. Need to be fixed. bK
344!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
345!> @todo chemistry error messages
346!
347!------------------------------------------------------------------------------!
348
349 MODULE chemistry_model_mod
350
351    USE advec_s_pw_mod,                                                                            &
352         ONLY:  advec_s_pw
353
354    USE advec_s_up_mod,                                                                            &
355         ONLY:  advec_s_up
356
357    USE advec_ws,                                                                                  &
358         ONLY:  advec_s_ws, ws_init_flags_scalar
359
360    USE diffusion_s_mod,                                                                           &
361         ONLY:  diffusion_s
362
363    USE kinds,                                                                                     &
364         ONLY:  iwp, wp
365
366    USE indices,                                                                                   &
367         ONLY:  advc_flags_s,                                                                      &
368                nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt, wall_flags_0
369
370    USE pegrid,                                                                                    &
371         ONLY: myid, threads_per_task
372
373    USE bulk_cloud_model_mod,                                                                      &
374         ONLY:  bulk_cloud_model
375
376    USE control_parameters,                                                                        &
377         ONLY:  bc_lr_cyc, bc_ns_cyc,                                                              &
378                bc_dirichlet_l,                                                                    &
379                bc_dirichlet_n,                                                                    &
380                bc_dirichlet_r,                                                                    &
381                bc_dirichlet_s,                                                                    &
382                bc_radiation_l,                                                                    &
383                bc_radiation_n,                                                                    &
384                bc_radiation_r,                                                                    &
385                bc_radiation_s,                                                                    &
386                debug_output,                                                                      &
387                dt_3d, humidity, initializing_actions, message_string,                             &
388                omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,          &
389                max_pr_user,                                                                       &
390                monotonic_limiter_z,                                                               &
391                scalar_advec,                                                                      &
392                timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
393
394    USE arrays_3d,                                                                                 &
395         ONLY:  exner, hyp, pt, q, ql, rdf_sc, tend, zu
396
397    USE chem_gasphase_mod,                                                                         &
398         ONLY:  atol, chem_gasphase_integrate, cs_mech, get_mechanism_name, nkppctrl,              &
399         nmaxfixsteps, nphot, nreact, nspec, nvar, phot_names, rtol, spc_names, t_steps, vl_dim
400
401    USE chem_modules
402
403    USE chem_photolysis_mod,                                                                       &
404        ONLY:  photolysis_control
405
406    USE cpulog,                                                                                    &
407        ONLY:  cpu_log, log_point_s
408
409    USE statistics
410
411    USE surface_mod,                                                                               &
412         ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
413
414    IMPLICIT NONE
415
416    PRIVATE
417    SAVE
418
419    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
420    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
421    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
422    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_av !< averaged concentrations of chemical species       
423    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
424                                                                            !< (only 1 timelevel required)
425    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
426                                                                            !< (e.g. solver type)
427    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
428                                                                            !< (e.g starting internal timestep of solver)
429!
430!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
431!-- approproate for chemicals compounds since they may accumulate too much.
432!-- If no proper boundary conditions from a DYNAMIC input file are available,
433!-- de-cycling applies the initial profiles at the inflow boundaries for
434!-- Dirichlet boundary conditions
435    LOGICAL ::  decycle_chem_lr    = .FALSE.    !< switch for setting decycling in left-right direction
436    LOGICAL ::  decycle_chem_ns    = .FALSE.    !< switch for setting decycling in south-norht direction
437    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
438         (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
439                              !< Decycling method at horizontal boundaries,
440                              !< 1=left, 2=right, 3=south, 4=north
441                              !< dirichlet = initial size distribution and
442                              !< chemical composition set for the ghost and
443                              !< first three layers
444                              !< neumann = zero gradient
445
446    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
447
448!
449!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
450    !
451    INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
452    INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
453    INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
454!
455!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
456    INTEGER(iwp) ::  ilu_grass              = 1       
457    INTEGER(iwp) ::  ilu_arable             = 2       
458    INTEGER(iwp) ::  ilu_permanent_crops    = 3         
459    INTEGER(iwp) ::  ilu_coniferous_forest  = 4         
460    INTEGER(iwp) ::  ilu_deciduous_forest   = 5         
461    INTEGER(iwp) ::  ilu_water_sea          = 6       
462    INTEGER(iwp) ::  ilu_urban              = 7       
463    INTEGER(iwp) ::  ilu_other              = 8 
464    INTEGER(iwp) ::  ilu_desert             = 9 
465    INTEGER(iwp) ::  ilu_ice                = 10 
466    INTEGER(iwp) ::  ilu_savanna            = 11 
467    INTEGER(iwp) ::  ilu_tropical_forest    = 12 
468    INTEGER(iwp) ::  ilu_water_inland       = 13 
469    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14 
470    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15 
471
472!
473!-- NH3/SO2 ratio regimes:
474    INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
475    INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
476    INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
477!
478!-- Default:
479    INTEGER, PARAMETER ::  iratns_default = iratns_low
480!
481!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
482    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
483         -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
484!
485!-- Set temperatures per land use for f_temp
486    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
487         12.0,  0.0, -999., 12.0,  8.0/)
488    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999.,      &
489         26.0, 20.0, -999., 26.0, 24.0 /)
490    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
491         40.0, 35.0, -999., 40.0, 39.0 /)
492!
493!-- Set f_min:
494    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01,  &
495         0.1, -999., 0.01, 0.04/)
496
497!
498!-- Set maximal conductance (m/s)
499!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
500    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
501         270., 150., -999., 300., 422./)/41000.
502!
503!-- Set max, min for vapour pressure deficit vpd
504    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3,      &
505         1.0, -999., 0.9, 2.8/) 
506    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999., 3.0, -999., -999., 3.0,     &
507         3.25, -999., 2.8, 4.5/) 
508
509    PUBLIC nreact
510    PUBLIC nspec               !< number of gas phase chemical species including constant compound (e.g. N2)
511    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
512    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
513    PUBLIC spec_conc_2 
514!   
515!-- Interface section
516    INTERFACE chem_3d_data_averaging
517       MODULE PROCEDURE chem_3d_data_averaging 
518    END INTERFACE chem_3d_data_averaging
519
520    INTERFACE chem_boundary_conds
521       MODULE PROCEDURE chem_boundary_conds
522       MODULE PROCEDURE chem_boundary_conds_decycle
523    END INTERFACE chem_boundary_conds
524
525    INTERFACE chem_check_data_output
526       MODULE PROCEDURE chem_check_data_output
527    END INTERFACE chem_check_data_output
528
529    INTERFACE chem_data_output_2d
530       MODULE PROCEDURE chem_data_output_2d
531    END INTERFACE chem_data_output_2d
532
533    INTERFACE chem_data_output_3d
534       MODULE PROCEDURE chem_data_output_3d
535    END INTERFACE chem_data_output_3d
536
537    INTERFACE chem_data_output_mask
538       MODULE PROCEDURE chem_data_output_mask
539    END INTERFACE chem_data_output_mask
540
541    INTERFACE chem_check_data_output_pr
542       MODULE PROCEDURE chem_check_data_output_pr
543    END INTERFACE chem_check_data_output_pr
544
545    INTERFACE chem_check_parameters
546       MODULE PROCEDURE chem_check_parameters
547    END INTERFACE chem_check_parameters
548
549    INTERFACE chem_define_netcdf_grid
550       MODULE PROCEDURE chem_define_netcdf_grid
551    END INTERFACE chem_define_netcdf_grid
552
553    INTERFACE chem_header
554       MODULE PROCEDURE chem_header
555    END INTERFACE chem_header
556
557    INTERFACE chem_init_arrays
558       MODULE PROCEDURE chem_init_arrays
559    END INTERFACE chem_init_arrays
560
561    INTERFACE chem_init
562       MODULE PROCEDURE chem_init
563    END INTERFACE chem_init
564
565    INTERFACE chem_init_profiles
566       MODULE PROCEDURE chem_init_profiles
567    END INTERFACE chem_init_profiles
568
569    INTERFACE chem_integrate
570       MODULE PROCEDURE chem_integrate_ij
571    END INTERFACE chem_integrate
572
573    INTERFACE chem_parin
574       MODULE PROCEDURE chem_parin
575    END INTERFACE chem_parin
576
577    INTERFACE chem_actions
578       MODULE PROCEDURE chem_actions
579       MODULE PROCEDURE chem_actions_ij
580    END INTERFACE chem_actions
581
582    INTERFACE chem_non_advective_processes
583       MODULE PROCEDURE chem_non_advective_processes
584       MODULE PROCEDURE chem_non_advective_processes_ij
585    END INTERFACE chem_non_advective_processes
586   
587    INTERFACE chem_exchange_horiz_bounds
588       MODULE PROCEDURE chem_exchange_horiz_bounds
589    END INTERFACE chem_exchange_horiz_bounds   
590
591    INTERFACE chem_prognostic_equations
592       MODULE PROCEDURE chem_prognostic_equations
593       MODULE PROCEDURE chem_prognostic_equations_ij
594    END INTERFACE chem_prognostic_equations
595
596    INTERFACE chem_rrd_local
597       MODULE PROCEDURE chem_rrd_local
598    END INTERFACE chem_rrd_local
599
600    INTERFACE chem_statistics
601       MODULE PROCEDURE chem_statistics
602    END INTERFACE chem_statistics
603
604    INTERFACE chem_swap_timelevel
605       MODULE PROCEDURE chem_swap_timelevel
606    END INTERFACE chem_swap_timelevel
607
608    INTERFACE chem_wrd_local
609       MODULE PROCEDURE chem_wrd_local 
610    END INTERFACE chem_wrd_local
611
612    INTERFACE chem_depo
613       MODULE PROCEDURE chem_depo 
614    END INTERFACE chem_depo
615
616    INTERFACE drydepos_gas_depac
617       MODULE PROCEDURE drydepos_gas_depac 
618    END INTERFACE drydepos_gas_depac
619
620    INTERFACE rc_special
621       MODULE PROCEDURE rc_special 
622    END INTERFACE rc_special
623
624    INTERFACE  rc_gw
625       MODULE PROCEDURE rc_gw 
626    END INTERFACE rc_gw
627
628    INTERFACE rw_so2 
629       MODULE PROCEDURE rw_so2 
630    END INTERFACE rw_so2
631
632    INTERFACE rw_nh3_sutton
633       MODULE PROCEDURE rw_nh3_sutton 
634    END INTERFACE rw_nh3_sutton
635
636    INTERFACE rw_constant
637       MODULE PROCEDURE rw_constant 
638    END INTERFACE rw_constant
639
640    INTERFACE rc_gstom
641       MODULE PROCEDURE rc_gstom 
642    END INTERFACE rc_gstom
643
644    INTERFACE rc_gstom_emb
645       MODULE PROCEDURE rc_gstom_emb 
646    END INTERFACE rc_gstom_emb
647
648    INTERFACE par_dir_diff
649       MODULE PROCEDURE par_dir_diff 
650    END INTERFACE par_dir_diff
651
652    INTERFACE rc_get_vpd
653       MODULE PROCEDURE rc_get_vpd 
654    END INTERFACE rc_get_vpd
655
656    INTERFACE rc_gsoil_eff
657       MODULE PROCEDURE rc_gsoil_eff 
658    END INTERFACE rc_gsoil_eff
659
660    INTERFACE rc_rinc
661       MODULE PROCEDURE rc_rinc 
662    END INTERFACE rc_rinc
663
664    INTERFACE rc_rctot
665       MODULE PROCEDURE rc_rctot 
666    END INTERFACE rc_rctot
667
668!    INTERFACE rc_comp_point_rc_eff
669!       MODULE PROCEDURE rc_comp_point_rc_eff
670!    END INTERFACE rc_comp_point_rc_eff
671
672    INTERFACE drydepo_aero_zhang_vd
673       MODULE PROCEDURE drydepo_aero_zhang_vd 
674    END INTERFACE drydepo_aero_zhang_vd
675
676    INTERFACE get_rb_cell
677       MODULE PROCEDURE  get_rb_cell
678    END INTERFACE get_rb_cell
679
680
681
682    PUBLIC chem_3d_data_averaging, chem_boundary_conds,                       &
683            chem_boundary_conds_decycle, chem_check_data_output,              &
684         chem_check_data_output_pr, chem_check_parameters,                    &
685         chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
686         chem_define_netcdf_grid, chem_header, chem_init, chem_init_arrays,   &
687         chem_init_profiles, chem_integrate, chem_parin,                      &
688         chem_actions, chem_prognostic_equations, chem_rrd_local,             &
689         chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo,     &
690         chem_non_advective_processes, chem_exchange_horiz_bounds
691
692 CONTAINS
693
694
695!------------------------------------------------------------------------------!
696! Description:
697! ------------
698!> Subroutine for averaging 3D data of chemical species. Due to the fact that
699!> the averaged chem arrays are allocated in chem_init, no if-query concerning
700!> the allocation is required (in any mode). Attention: If you just specify an
701!> averaged output quantity in the _p3dr file during restarts the first output
702!> includes the time between the beginning of the restart run and the first
703!> output time (not necessarily the whole averaging_interval you have
704!> specified in your _p3d/_p3dr file )
705!------------------------------------------------------------------------------!
706 SUBROUTINE chem_3d_data_averaging( mode, variable )
707
708
709    USE control_parameters
710
711    CHARACTER (LEN=*) ::  mode     !<
712    CHARACTER (LEN=*) ::  variable !<
713
714    LOGICAL ::  match_def  !< flag indicating default-type surface
715    LOGICAL ::  match_lsm  !< flag indicating natural-type surface
716    LOGICAL ::  match_usm  !< flag indicating urban-type surface
717
718    INTEGER(iwp) ::  i                  !< grid index x direction
719    INTEGER(iwp) ::  j                  !< grid index y direction
720    INTEGER(iwp) ::  k                  !< grid index z direction
721    INTEGER(iwp) ::  m                  !< running index surface type
722    INTEGER(iwp) ::  lsp               !< running index for chem spcs
723
724    IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_')  )  THEN
725
726       IF ( mode == 'allocate' )  THEN
727
728          DO  lsp = 1, nspec
729             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
730                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
731                chem_species(lsp)%conc_av = 0.0_wp
732             ENDIF
733          ENDDO
734
735       ELSEIF ( mode == 'sum' )  THEN
736
737          DO  lsp = 1, nspec
738             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
739                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
740                DO  i = nxlg, nxrg
741                   DO  j = nysg, nyng
742                      DO  k = nzb, nzt+1
743                         chem_species(lsp)%conc_av(k,j,i) =                              &
744                                           chem_species(lsp)%conc_av(k,j,i) +            &
745                                           chem_species(lsp)%conc(k,j,i)
746                      ENDDO
747                   ENDDO
748                ENDDO
749             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
750                DO  i = nxl, nxr
751                   DO  j = nys, nyn
752                      match_def = surf_def_h(0)%start_index(j,i) <=                      &
753                           surf_def_h(0)%end_index(j,i)
754                      match_lsm = surf_lsm_h%start_index(j,i) <=                         &
755                           surf_lsm_h%end_index(j,i)
756                      match_usm = surf_usm_h%start_index(j,i) <=                         &
757                           surf_usm_h%end_index(j,i)
758
759                      IF ( match_def )  THEN
760                         m = surf_def_h(0)%end_index(j,i)
761                         chem_species(lsp)%cssws_av(j,i) =                               &
762                              chem_species(lsp)%cssws_av(j,i) + &
763                              surf_def_h(0)%cssws(lsp,m)
764                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
765                         m = surf_lsm_h%end_index(j,i)
766                         chem_species(lsp)%cssws_av(j,i) =                               &
767                              chem_species(lsp)%cssws_av(j,i) + &
768                              surf_lsm_h%cssws(lsp,m)
769                      ELSEIF ( match_usm )  THEN
770                         m = surf_usm_h%end_index(j,i)
771                         chem_species(lsp)%cssws_av(j,i) =                               &
772                              chem_species(lsp)%cssws_av(j,i) + &
773                              surf_usm_h%cssws(lsp,m)
774                      ENDIF
775                   ENDDO
776                ENDDO
777             ENDIF
778          ENDDO
779
780       ELSEIF ( mode == 'average' )  THEN
781
782          DO  lsp = 1, nspec
783             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
784                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
785                DO  i = nxlg, nxrg
786                   DO  j = nysg, nyng
787                      DO  k = nzb, nzt+1
788                         chem_species(lsp)%conc_av(k,j,i) =                              &
789                             chem_species(lsp)%conc_av(k,j,i) /                          &
790                             REAL( average_count_3d, KIND=wp )
791                      ENDDO
792                   ENDDO
793                ENDDO
794
795             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
796                DO  i = nxlg, nxrg
797                   DO  j = nysg, nyng
798                      chem_species(lsp)%cssws_av(j,i) =                                  &
799                      chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
800                   ENDDO
801                ENDDO
802                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
803             ENDIF
804          ENDDO
805       ENDIF
806
807    ENDIF
808
809 END SUBROUTINE chem_3d_data_averaging
810
811   
812!------------------------------------------------------------------------------!
813! Description:
814! ------------
815!> Subroutine to initialize and set all boundary conditions for chemical species
816!------------------------------------------------------------------------------!
817 SUBROUTINE chem_boundary_conds( mode )                                           
818
819    USE control_parameters,                                                    & 
820        ONLY:  bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s
821
822    USE arrays_3d,                                                             &     
823        ONLY:  dzu                                               
824
825    USE surface_mod,                                                           &
826        ONLY:  bc_h                                                           
827
828    CHARACTER (LEN=*), INTENT(IN) ::  mode
829    INTEGER(iwp) ::  i                            !< grid index x direction.
830    INTEGER(iwp) ::  j                            !< grid index y direction.
831    INTEGER(iwp) ::  k                            !< grid index z direction.
832    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
833    INTEGER(iwp) ::  m                            !< running index surface elements.
834    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
835
836
837    SELECT CASE ( TRIM( mode ) )       
838       CASE ( 'init' )
839
840          IF ( bc_cs_b == 'dirichlet' )  THEN
841             ibc_cs_b = 0 
842          ELSEIF ( bc_cs_b == 'neumann' )  THEN
843             ibc_cs_b = 1 
844          ELSE
845             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
846             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
847          ENDIF                                                                 
848!
849!--       Set Integer flags and check for possible erroneous settings for top
850!--       boundary condition.
851          IF ( bc_cs_t == 'dirichlet' )  THEN
852             ibc_cs_t = 0 
853          ELSEIF ( bc_cs_t == 'neumann' )  THEN
854             ibc_cs_t = 1
855          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
856             ibc_cs_t = 2
857          ELSEIF ( bc_cs_t == 'nested' )  THEN         
858             ibc_cs_t = 3
859          ELSE
860             message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'     
861             CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
862          ENDIF
863
864       CASE ( 'set_bc_bottomtop' )                   
865!
866!--       Boundary condtions for chemical species at horizontal walls     
867          DO  lsp = 1, nspec                                                     
868             IF ( ibc_cs_b == 0 )  THEN
869                DO  l = 0, 1 
870                   !$OMP PARALLEL DO PRIVATE( i, j, k )
871                   DO  m = 1, bc_h(l)%ns
872                       i = bc_h(l)%i(m)           
873                       j = bc_h(l)%j(m)
874                       k = bc_h(l)%k(m)
875                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
876                                      chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 
877                   ENDDO                                       
878                ENDDO                                       
879
880             ELSEIF ( ibc_cs_b == 1 )  THEN
881!
882!--             in boundary_conds there is som extra loop over m here for passive tracer
883                DO  l = 0, 1
884                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
885                   DO m = 1, bc_h(l)%ns
886                      i = bc_h(l)%i(m)           
887                      j = bc_h(l)%j(m)
888                      k = bc_h(l)%k(m)
889                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
890                                         chem_species(lsp)%conc_p(k,j,i)
891
892                   ENDDO
893                ENDDO
894             ENDIF
895       ENDDO    ! end lsp loop 
896!
897!--    Top boundary conditions for chemical species - Should this not be done for all species?
898          IF ( ibc_cs_t == 0 )  THEN                   
899             DO  lsp = 1, nspec
900                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
901             ENDDO
902          ELSEIF ( ibc_cs_t == 1 )  THEN
903             DO  lsp = 1, nspec
904                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
905             ENDDO
906          ELSEIF ( ibc_cs_t == 2 )  THEN
907             DO  lsp = 1, nspec
908                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
909             ENDDO
910          ENDIF
911
912       CASE ( 'set_bc_lateral' )                       
913!
914!--             Lateral boundary conditions for chem species at inflow boundary
915!--             are automatically set when chem_species concentration is
916!--             initialized. The initially set value at the inflow boundary is not
917!--             touched during time integration, hence, this boundary value remains
918!--             at a constant value, which is the concentration that flows into the
919!--             domain.                                                           
920!--             Lateral boundary conditions for chem species at outflow boundary
921
922          IF ( bc_radiation_s )  THEN
923             DO  lsp = 1, nspec
924                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
925             ENDDO
926          ELSEIF ( bc_radiation_n )  THEN
927             DO  lsp = 1, nspec
928                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
929             ENDDO
930          ELSEIF ( bc_radiation_l )  THEN
931             DO  lsp = 1, nspec
932                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
933             ENDDO
934          ELSEIF ( bc_radiation_r )  THEN
935             DO  lsp = 1, nspec
936                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
937             ENDDO
938          ENDIF
939
940    END SELECT
941
942 END SUBROUTINE chem_boundary_conds
943
944
945!------------------------------------------------------------------------------!
946! Description:
947! ------------
948!> Boundary conditions for prognostic variables in chemistry: decycling in the
949!> x-direction-
950!> Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
951!> approproate for chemicals compounds since they may accumulate too much.
952!> If no proper boundary conditions from a DYNAMIC input file are available,
953!> de-cycling applies the initial profiles at the inflow boundaries for
954!> Dirichlet boundary conditions
955!------------------------------------------------------------------------------!
956 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
957
958
959   INTEGER(iwp) ::  boundary  !<
960   INTEGER(iwp) ::  ee        !<
961   INTEGER(iwp) ::  copied    !<
962   INTEGER(iwp) ::  i         !<
963   INTEGER(iwp) ::  j         !<
964   INTEGER(iwp) ::  k         !<
965   INTEGER(iwp) ::  ss        !<
966
967   REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
968   REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
969   REAL(wp) ::  flag          !< flag to mask topography grid points
970
971
972   flag = 0.0_wp
973   !
974   !-- Left and right boundaries
975   IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
976
977      DO  boundary = 1, 2
978
979         IF ( decycle_method(boundary) == 'dirichlet' )  THEN
980            !
981            !--          Initial profile is copied to ghost and first three layers
982            ss = 1
983            ee = 0
984            IF ( boundary == 1  .AND.  nxl == 0 )  THEN
985               ss = nxlg
986               ee = nxl-1
987            ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
988               ss = nxr+1
989               ee = nxrg
990            ENDIF
991
992            DO  i = ss, ee
993               DO  j = nysg, nyng
994                  DO  k = nzb+1, nzt
995                     flag = MERGE( 1.0_wp, 0.0_wp,                            &
996                          BTEST( wall_flags_0(k,j,i), 0 ) )
997                     cs_3d(k,j,i) = cs_pr_init(k) * flag
998                  ENDDO
999               ENDDO
1000            ENDDO
1001
1002         ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
1003            !
1004            !--          The value at the boundary is copied to the ghost layers to simulate
1005            !--          an outlet with zero gradient
1006            ss = 1
1007            ee = 0
1008            IF ( boundary == 1  .AND.  nxl == 0 )  THEN
1009               ss = nxlg
1010               ee = nxl-1
1011               copied = nxl
1012            ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
1013               ss = nxr+1
1014               ee = nxrg
1015               copied = nxr
1016            ENDIF
1017
1018            DO  i = ss, ee
1019               DO  j = nysg, nyng
1020                  DO  k = nzb+1, nzt
1021                     flag = MERGE( 1.0_wp, 0.0_wp,                            &
1022                          BTEST( wall_flags_0(k,j,i), 0 ) )
1023                     cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
1024                  ENDDO
1025               ENDDO
1026            ENDDO
1027
1028         ELSE
1029            WRITE(message_string,*)                                           &
1030                 'unknown decycling method: decycle_method (', &
1031                 boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
1032            CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
1033                 1, 2, 0, 6, 0 )
1034         ENDIF
1035      ENDDO
1036   ENDIF
1037   !
1038   !-- South and north boundaries
1039   IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
1040
1041      DO  boundary = 3, 4
1042
1043         IF ( decycle_method(boundary) == 'dirichlet' )  THEN
1044            !
1045            !--          Initial profile is copied to ghost and first three layers
1046            ss = 1
1047            ee = 0
1048            IF ( boundary == 3  .AND.  nys == 0 )  THEN
1049               ss = nysg
1050               ee = nys-1
1051            ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
1052               ss = nyn+1
1053               ee = nyng
1054            ENDIF
1055
1056            DO  i = nxlg, nxrg
1057               DO  j = ss, ee
1058                  DO  k = nzb+1, nzt
1059                     flag = MERGE( 1.0_wp, 0.0_wp,                            &
1060                          BTEST( wall_flags_0(k,j,i), 0 ) )
1061                     cs_3d(k,j,i) = cs_pr_init(k) * flag
1062                  ENDDO
1063               ENDDO
1064            ENDDO
1065
1066
1067         ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
1068            !
1069            !--          The value at the boundary is copied to the ghost layers to simulate
1070            !--          an outlet with zero gradient
1071            ss = 1
1072            ee = 0
1073            IF ( boundary == 3  .AND.  nys == 0 )  THEN
1074               ss = nysg
1075               ee = nys-1
1076               copied = nys
1077            ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
1078               ss = nyn+1
1079               ee = nyng
1080               copied = nyn
1081            ENDIF
1082
1083            DO  i = nxlg, nxrg
1084               DO  j = ss, ee
1085                  DO  k = nzb+1, nzt
1086                     flag = MERGE( 1.0_wp, 0.0_wp,                            &
1087                          BTEST( wall_flags_0(k,j,i), 0 ) )
1088                     cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
1089                  ENDDO
1090               ENDDO
1091            ENDDO
1092
1093         ELSE
1094            WRITE(message_string,*)                                           &
1095                 'unknown decycling method: decycle_method (', &
1096                 boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
1097            CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
1098                 1, 2, 0, 6, 0 )
1099         ENDIF
1100      ENDDO
1101   ENDIF
1102
1103
1104 END SUBROUTINE chem_boundary_conds_decycle
1105
1106
1107!------------------------------------------------------------------------------!
1108! Description:
1109! ------------
1110!> Subroutine for checking data output for chemical species
1111!------------------------------------------------------------------------------!
1112 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
1113
1114
1115    CHARACTER (LEN=*) ::  unit     !<
1116    CHARACTER (LEN=*) ::  var      !<
1117
1118    INTEGER(iwp) ::  i
1119    INTEGER(iwp) ::  lsp
1120    INTEGER(iwp) ::  ilen
1121    INTEGER(iwp) ::  k
1122
1123    CHARACTER(LEN=16)    ::  spec_name
1124
1125!
1126!-- Next statement is to avoid compiler warnings about unused variables
1127    IF ( ( i + ilen + k ) > 0  .OR.  var(1:1) == ' ' )  CONTINUE
1128
1129    unit = 'illegal'
1130
1131    spec_name = TRIM( var(4:) )             !< var 1:3 is 'kc_' or 'em_'.
1132
1133    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
1134       DO  lsp=1,nspec
1135          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1136             unit = 'mol m-2 s-1'
1137          ENDIF
1138!
1139!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code
1140!--       as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1141!--       set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
1142          IF (spec_name(1:2) == 'PM')  THEN
1143             unit = 'kg m-2 s-1'
1144          ENDIF
1145       ENDDO
1146
1147    ELSE
1148
1149       DO  lsp=1,nspec
1150          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1151             unit = 'ppm'
1152          ENDIF
1153!
1154!--            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
1155!--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1156!--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
1157          IF (spec_name(1:2) == 'PM')  THEN 
1158            unit = 'kg m-3'
1159          ENDIF
1160       ENDDO
1161
1162       DO  lsp=1,nphot
1163          IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
1164             unit = 'sec-1'
1165          ENDIF
1166       ENDDO
1167    ENDIF
1168
1169
1170    RETURN
1171 END SUBROUTINE chem_check_data_output
1172
1173
1174!------------------------------------------------------------------------------!
1175! Description:
1176! ------------
1177!> Subroutine for checking data output of profiles for chemistry model
1178!------------------------------------------------------------------------------!
1179 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1180
1181    USE arrays_3d
1182
1183    USE control_parameters,                                                    &
1184        ONLY:  data_output_pr, message_string
1185
1186    USE profil_parameter
1187
1188    USE statistics
1189
1190
1191    CHARACTER (LEN=*) ::  unit     !<
1192    CHARACTER (LEN=*) ::  variable !<
1193    CHARACTER (LEN=*) ::  dopr_unit
1194    CHARACTER (LEN=16) ::  spec_name
1195
1196    INTEGER(iwp) ::  var_count, lsp  !<
1197
1198
1199    spec_name = TRIM( variable(4:) )   
1200
1201       IF (  .NOT.  air_chemistry )  THEN
1202          message_string = 'data_output_pr = ' //                        &
1203          TRIM( data_output_pr(var_count) ) // ' is not imp' // &
1204                      'lemented for air_chemistry = .FALSE.'
1205          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
1206
1207       ELSE
1208          DO  lsp = 1, nspec
1209             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1210                cs_pr_count = cs_pr_count+1
1211                cs_pr_index(cs_pr_count) = lsp
1212                dopr_index(var_count) = pr_palm+cs_pr_count 
1213                dopr_unit  = 'ppm'
1214                IF (spec_name(1:2) == 'PM')  THEN
1215                     dopr_unit = 'kg m-3'
1216                ENDIF
1217                   hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
1218                   unit = dopr_unit 
1219             ENDIF
1220          ENDDO
1221       ENDIF
1222
1223 END SUBROUTINE chem_check_data_output_pr
1224
1225
1226!------------------------------------------------------------------------------!
1227! Description:
1228! ------------
1229!> Check parameters routine for chemistry_model_mod
1230!------------------------------------------------------------------------------!
1231 SUBROUTINE chem_check_parameters
1232
1233
1234    LOGICAL  ::  found
1235    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
1236    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
1237!
1238!-- check for chemical reactions status
1239    IF ( chem_gasphase_on )  THEN
1240       message_string = 'Chemical reactions: ON'
1241       CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
1242    ELSEIF ( .NOT. (chem_gasphase_on) )  THEN
1243       message_string = 'Chemical reactions: OFF'
1244       CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
1245    ENDIF
1246!
1247!-- check for chemistry time-step
1248    IF ( call_chem_at_all_substeps )  THEN
1249       message_string = 'Chemistry is calculated at all meteorology time-step'
1250       CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
1251    ELSEIF ( .not. (call_chem_at_all_substeps) )  THEN
1252       message_string = 'Sub-time-steps are skipped for chemistry time-steps'
1253       CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
1254    ENDIF
1255!
1256!-- check for photolysis scheme
1257    IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
1258       message_string = 'Incorrect photolysis scheme selected, please check spelling'
1259       CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
1260    ENDIF
1261!
1262!-- check for  decycling of chem species
1263    IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) )  THEN
1264       message_string = 'Incorrect boundary conditions. Only neumann or '   &
1265                // 'dirichlet &available for decycling chemical species '
1266       CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
1267    ENDIF
1268!
1269!-- check for chemical mechanism used
1270    CALL get_mechanism_name
1271    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
1272       message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
1273       CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
1274    ENDIF
1275!
1276!-- chem_check_parameters is called before the array chem_species is allocated!
1277!-- temporary switch of this part of the check
1278!    RETURN                !bK commented
1279
1280    CALL chem_init_internal
1281!
1282!-- check for initial chem species input
1283    lsp_usr = 1
1284    lsp     = 1
1285    DO WHILE ( cs_name (lsp_usr) /= 'novalue')
1286       found = .FALSE.
1287       DO  lsp = 1, nvar
1288          IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) )  THEN
1289             found = .TRUE.
1290             EXIT
1291          ENDIF
1292       ENDDO
1293       IF ( .NOT.  found )  THEN
1294          message_string = 'Unused/incorrect input for initial surface value: ' //     &
1295                            TRIM( cs_name(lsp_usr) )
1296          CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 )
1297       ENDIF
1298       lsp_usr = lsp_usr + 1
1299    ENDDO
1300!
1301!-- check for surface  emission flux chem species
1302    lsp_usr = 1
1303    lsp     = 1
1304    DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
1305       found = .FALSE.
1306       DO  lsp = 1, nvar
1307          IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) )  THEN
1308             found = .TRUE.
1309             EXIT
1310          ENDIF
1311       ENDDO
1312       IF ( .NOT.  found )  THEN
1313          message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
1314                            // TRIM( surface_csflux_name(lsp_usr) )
1315          CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 )
1316       ENDIF
1317       lsp_usr = lsp_usr + 1
1318    ENDDO
1319
1320 END SUBROUTINE chem_check_parameters
1321
1322
1323!------------------------------------------------------------------------------!
1324! Description:
1325! ------------
1326!> Subroutine defining 2D output variables for chemical species
1327!> @todo: Remove "mode" from argument list, not used.
1328!------------------------------------------------------------------------------!
1329 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
1330                                  two_d, nzb_do, nzt_do, fill_value )
1331
1332
1333    CHARACTER (LEN=*) ::  grid       !<
1334    CHARACTER (LEN=*) ::  mode       !<
1335    CHARACTER (LEN=*) ::  variable   !<
1336    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
1337    INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
1338    INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
1339    LOGICAL      ::  found           !<
1340    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
1341    REAL(wp)     ::  fill_value 
1342    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf 
1343
1344!
1345!-- local variables.
1346    CHARACTER(LEN=16)    ::  spec_name
1347    INTEGER(iwp) ::  lsp
1348    INTEGER(iwp) ::  i               !< grid index along x-direction
1349    INTEGER(iwp) ::  j               !< grid index along y-direction
1350    INTEGER(iwp) ::  k               !< grid index along z-direction
1351    INTEGER(iwp) ::  m               !< running indices for surfaces
1352    INTEGER(iwp) ::  char_len        !< length of a character string
1353!
1354!-- Next statement is to avoid compiler warnings about unused variables
1355    IF ( mode(1:1) == ' '  .OR.  two_d )  CONTINUE
1356
1357    found = .FALSE.
1358    char_len  = LEN_TRIM( variable )
1359
1360    spec_name = TRIM( variable(4:char_len-3) )
1361!
1362!-- Output of emission values, i.e. surface fluxes cssws.
1363    IF ( variable(1:3) == 'em_' )  THEN
1364
1365       local_pf = 0.0_wp
1366
1367       DO  lsp = 1, nvar
1368          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1369!
1370!--          No average output for now.
1371             DO  m = 1, surf_lsm_h%ns
1372                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =              &
1373                              local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1)  &
1374                                            + surf_lsm_h%cssws(lsp,m)
1375             ENDDO
1376             DO  m = 1, surf_usm_h%ns
1377                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =           &
1378                              local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1)  &
1379                                            + surf_usm_h%cssws(lsp,m)
1380             ENDDO
1381             grid = 'zu'
1382             found = .TRUE.
1383          ENDIF
1384       ENDDO
1385
1386    ELSE
1387
1388       DO  lsp=1,nspec
1389          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                           &
1390                ( (variable(char_len-2:) == '_xy')  .OR.                                           &
1391                  (variable(char_len-2:) == '_xz')  .OR.                                           &
1392                  (variable(char_len-2:) == '_yz') ) )  THEN             
1393!
1394!--   todo: remove or replace by "CALL message" mechanism (kanani)
1395!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
1396!                                                             TRIM( chem_species(lsp)%name )       
1397             IF (av == 0)  THEN
1398                DO  i = nxl, nxr
1399                   DO  j = nys, nyn
1400                      DO  k = nzb_do, nzt_do
1401                           local_pf(i,j,k) = MERGE(                                                &
1402                                              chem_species(lsp)%conc(k,j,i),                       &
1403                                              REAL( fill_value, KIND = wp ),                       &
1404                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1405                      ENDDO
1406                   ENDDO
1407                ENDDO
1408
1409             ELSE
1410                DO  i = nxl, nxr
1411                   DO  j = nys, nyn
1412                      DO  k = nzb_do, nzt_do
1413                           local_pf(i,j,k) = MERGE(                                                &
1414                                              chem_species(lsp)%conc_av(k,j,i),                    &
1415                                              REAL( fill_value, KIND = wp ),                       &
1416                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1417                      ENDDO
1418                   ENDDO
1419                ENDDO
1420             ENDIF
1421             grid = 'zu'           
1422             found = .TRUE.
1423          ENDIF
1424       ENDDO
1425    ENDIF
1426
1427    RETURN
1428
1429 END SUBROUTINE chem_data_output_2d     
1430
1431
1432!------------------------------------------------------------------------------!
1433! Description:
1434! ------------
1435!> Subroutine defining 3D output variables for chemical species
1436!------------------------------------------------------------------------------!
1437 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1438
1439
1440    USE surface_mod
1441
1442    CHARACTER (LEN=*)    ::  variable     !<
1443    INTEGER(iwp)         ::  av           !<
1444    INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
1445    INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
1446
1447    LOGICAL      ::  found                !<
1448
1449    REAL(wp)             ::  fill_value   !<
1450    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
1451!
1452!-- local variables
1453    CHARACTER(LEN=16)    ::  spec_name
1454    INTEGER(iwp)         ::  i
1455    INTEGER(iwp)         ::  j
1456    INTEGER(iwp)         ::  k
1457    INTEGER(iwp)         ::  m       !< running indices for surfaces
1458    INTEGER(iwp)         ::  l
1459    INTEGER(iwp)         ::  lsp     !< running index for chem spcs
1460
1461
1462    found = .FALSE.
1463    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) )  THEN
1464       RETURN
1465    ENDIF
1466
1467    spec_name = TRIM( variable(4:) )
1468
1469    IF ( variable(1:3) == 'em_' )  THEN
1470
1471       DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
1472          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1473         
1474             local_pf = 0.0_wp
1475!
1476!--          no average for now
1477             DO  m = 1, surf_usm_h%ns
1478                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
1479                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
1480             ENDDO
1481             DO  m = 1, surf_lsm_h%ns
1482                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
1483                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
1484             ENDDO
1485             DO  l = 0, 3
1486                DO  m = 1, surf_usm_v(l)%ns
1487                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
1488                     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)
1489                ENDDO
1490                DO  m = 1, surf_lsm_v(l)%ns
1491                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
1492                      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)
1493                ENDDO
1494             ENDDO
1495             found = .TRUE.
1496          ENDIF
1497       ENDDO
1498    ELSE
1499      DO  lsp = 1, nspec
1500         IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1501!
1502!--   todo: remove or replace by "CALL message" mechanism (kanani)
1503!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1504!                                                           TRIM( chem_species(lsp)%name )       
1505            IF (av == 0)  THEN
1506               DO  i = nxl, nxr
1507                  DO  j = nys, nyn
1508                     DO  k = nzb_do, nzt_do
1509                         local_pf(i,j,k) = MERGE(                             &
1510                                             chem_species(lsp)%conc(k,j,i),   &
1511                                             REAL( fill_value, KIND = wp ),   &
1512                                             BTEST( wall_flags_0(k,j,i), 0 ) )
1513                     ENDDO
1514                  ENDDO
1515               ENDDO
1516
1517            ELSE
1518
1519               DO  i = nxl, nxr
1520                  DO  j = nys, nyn
1521                     DO  k = nzb_do, nzt_do
1522                         local_pf(i,j,k) = MERGE(                             &
1523                                             chem_species(lsp)%conc_av(k,j,i),&
1524                                             REAL( fill_value, KIND = wp ),   &
1525                                             BTEST( wall_flags_0(k,j,i), 0 ) )
1526                     ENDDO
1527                  ENDDO
1528               ENDDO
1529            ENDIF
1530            found = .TRUE.
1531         ENDIF
1532      ENDDO
1533    ENDIF
1534
1535    RETURN
1536
1537 END SUBROUTINE chem_data_output_3d
1538
1539
1540!------------------------------------------------------------------------------!
1541! Description:
1542! ------------
1543!> Subroutine defining mask output variables for chemical species
1544!------------------------------------------------------------------------------!
1545 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
1546
1547
1548    USE control_parameters
1549
1550    USE surface_mod,                                                                  &
1551        ONLY:  get_topography_top_index_ji
1552
1553
1554    CHARACTER(LEN=5)  ::  grid        !< flag to distinquish between staggered grids
1555    CHARACTER(LEN=16) ::  spec_name
1556    CHARACTER(LEN=*)  ::  variable    !<
1557   
1558    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
1559    INTEGER(iwp) ::  lsp
1560    INTEGER(iwp) ::  i               !< grid index along x-direction
1561    INTEGER(iwp) ::  j               !< grid index along y-direction
1562    INTEGER(iwp) ::  k               !< grid index along z-direction
1563    INTEGER(iwp) ::  mid             !< masked output running index
1564    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
1565   
1566    LOGICAL ::  found
1567   
1568    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
1569              local_pf   !<
1570!
1571!-- local variables.
1572
1573    found = .TRUE.
1574    grid = 's'
1575
1576    spec_name = TRIM( variable(4:) )
1577
1578    DO  lsp=1,nspec
1579       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN             
1580!
1581!-- todo: remove or replace by "CALL message" mechanism (kanani)
1582!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1583!                                                        TRIM( chem_species(lsp)%name )       
1584          IF (av == 0)  THEN
1585             IF ( .NOT. mask_surface(mid) )  THEN
1586
1587                DO  i = 1, mask_size_l(mid,1)
1588                   DO  j = 1, mask_size_l(mid,2) 
1589                      DO  k = 1, mask_size(mid,3) 
1590                          local_pf(i,j,k) = chem_species(lsp)%conc(  &
1591                                               mask_k(mid,k),        &
1592                                               mask_j(mid,j),        &
1593                                               mask_i(mid,i)      )
1594                      ENDDO
1595                   ENDDO
1596                ENDDO
1597
1598             ELSE
1599!             
1600!--             Terrain-following masked output
1601                DO  i = 1, mask_size_l(mid,1)
1602                   DO  j = 1, mask_size_l(mid,2)
1603!             
1604!--                   Get k index of highest horizontal surface
1605                      topo_top_ind = get_topography_top_index_ji( &
1606                                        mask_j(mid,j),  &
1607                                        mask_i(mid,i),  &
1608                                        grid                    )
1609!             
1610!--                   Save output array
1611                      DO  k = 1, mask_size_l(mid,3)
1612                         local_pf(i,j,k) = chem_species(lsp)%conc( &
1613                                              MIN( topo_top_ind+mask_k(mid,k), &
1614                                                   nzt+1 ),        &
1615                                              mask_j(mid,j),       &
1616                                              mask_i(mid,i)      )
1617                      ENDDO
1618                   ENDDO
1619                ENDDO
1620
1621             ENDIF
1622          ELSE
1623             IF ( .NOT. mask_surface(mid) )  THEN
1624
1625                DO  i = 1, mask_size_l(mid,1)
1626                   DO  j = 1, mask_size_l(mid,2)
1627                      DO  k =  1, mask_size_l(mid,3)
1628                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
1629                                               mask_k(mid,k),           &
1630                                               mask_j(mid,j),           &
1631                                               mask_i(mid,i)         )
1632                      ENDDO
1633                   ENDDO
1634                ENDDO
1635
1636             ELSE
1637!             
1638!--             Terrain-following masked output
1639                DO  i = 1, mask_size_l(mid,1)
1640                   DO  j = 1, mask_size_l(mid,2)
1641!             
1642!--                   Get k index of highest horizontal surface
1643                      topo_top_ind = get_topography_top_index_ji( &
1644                                        mask_j(mid,j),  &
1645                                        mask_i(mid,i),  &
1646                                        grid                    )
1647!             
1648!--                   Save output array
1649                      DO  k = 1, mask_size_l(mid,3)
1650                         local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
1651                                              MIN( topo_top_ind+mask_k(mid,k), &
1652                                                   nzt+1 ),            &
1653                                              mask_j(mid,j),           &
1654                                              mask_i(mid,i)         )
1655                      ENDDO
1656                   ENDDO
1657                ENDDO
1658
1659             ENDIF
1660
1661          ENDIF
1662          found = .FALSE.
1663       ENDIF
1664    ENDDO
1665
1666    RETURN
1667
1668 END SUBROUTINE chem_data_output_mask     
1669
1670
1671!------------------------------------------------------------------------------!
1672! Description:
1673! ------------
1674!> Subroutine defining appropriate grid for netcdf variables.
1675!> It is called out from subroutine netcdf.
1676!------------------------------------------------------------------------------!
1677 SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1678
1679
1680    CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
1681    LOGICAL, INTENT(OUT)           ::  found        !<
1682    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
1683    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
1684    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
1685
1686    found  = .TRUE.
1687
1688    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                   !< always the same grid for chemistry variables
1689       grid_x = 'x'
1690       grid_y = 'y'
1691       grid_z = 'zu'                             
1692    ELSE
1693       found  = .FALSE.
1694       grid_x = 'none'
1695       grid_y = 'none'
1696       grid_z = 'none'
1697    ENDIF
1698
1699
1700 END SUBROUTINE chem_define_netcdf_grid
1701
1702
1703!------------------------------------------------------------------------------!
1704! Description:
1705! ------------
1706!> Subroutine defining header output for chemistry model
1707!------------------------------------------------------------------------------!
1708 SUBROUTINE chem_header( io )
1709
1710
1711    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1712    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
1713    INTEGER(iwp)  :: cs_fixed 
1714    CHARACTER (LEN=80)  :: docsflux_chr
1715    CHARACTER (LEN=80)  :: docsinit_chr
1716!
1717! Get name of chemical mechanism from chem_gasphase_mod
1718    CALL get_mechanism_name
1719!
1720!-- Write chemistry model  header
1721    WRITE( io, 1 )
1722!
1723!-- Gasphase reaction status
1724    IF ( chem_gasphase_on )  THEN
1725       WRITE( io, 2 )
1726    ELSE
1727       WRITE( io, 3 )
1728    ENDIF
1729!
1730!-- Chemistry time-step
1731    WRITE ( io, 4 ) cs_time_step
1732!
1733!-- Emission mode info
1734!-- At the moment the evaluation is done with both emiss_lod and mode_emis
1735!-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
1736!-- conditions can be removed (ecc 20190513)
1737    IF     ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') )        THEN
1738        WRITE ( io, 5 )
1739    ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') )  THEN
1740        WRITE ( io, 6 )
1741    ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') )  THEN
1742        WRITE ( io, 7 )
1743    ENDIF
1744!
1745!-- Photolysis scheme info
1746    IF ( photolysis_scheme == "simple" )  THEN
1747       WRITE( io, 8 ) 
1748    ELSEIF (photolysis_scheme == "constant" )  THEN
1749       WRITE( io, 9 )
1750    ENDIF
1751!
1752!-- Emission flux info
1753    lsp = 1
1754    docsflux_chr ='Chemical species for surface emission flux: ' 
1755    DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1756       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
1757       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
1758          WRITE ( io, 10 )  docsflux_chr
1759          docsflux_chr = '       '
1760       ENDIF
1761       lsp = lsp + 1
1762    ENDDO
1763
1764    IF ( docsflux_chr /= '' )  THEN
1765       WRITE ( io, 10 )  docsflux_chr
1766    ENDIF
1767!
1768!-- initializatoin of Surface and profile chemical species
1769
1770    lsp = 1
1771    docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
1772    DO WHILE ( cs_name(lsp) /= 'novalue' )
1773       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
1774       IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
1775        WRITE ( io, 11 )  docsinit_chr
1776        docsinit_chr = '       '
1777       ENDIF
1778       lsp = lsp + 1
1779    ENDDO
1780
1781    IF ( docsinit_chr /= '' )  THEN
1782       WRITE ( io, 11 )  docsinit_chr
1783    ENDIF
1784!
1785!-- number of variable and fix chemical species and number of reactions
1786    cs_fixed = nspec - nvar
1787
1788    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech 
1789    WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
1790    WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
1791    WRITE ( io, * ) '   --> Total number of reactions : ', nreact
1792
1793
17941   FORMAT (//' Chemistry model information:'/                                  &
1795           ' ----------------------------'/)
17962   FORMAT ('    --> Chemical reactions are turned on')
17973   FORMAT ('    --> Chemical reactions are turned off')
17984   FORMAT ('    --> Time-step for chemical species: ',F6.2, ' s')
17995   FORMAT ('    --> Emission mode = DEFAULT ')
18006   FORMAT ('    --> Emission mode = PARAMETERIZED ')
18017   FORMAT ('    --> Emission mode = PRE-PROCESSED ')
18028   FORMAT ('    --> Photolysis scheme used =  simple ')
18039   FORMAT ('    --> Photolysis scheme used =  constant ')
180410  FORMAT (/'    ',A) 
180511  FORMAT (/'    ',A) 
1806!
1807!
1808 END SUBROUTINE chem_header
1809
1810
1811!------------------------------------------------------------------------------!
1812! Description:
1813! ------------
1814!> Subroutine initializating chemistry_model_mod specific arrays
1815!------------------------------------------------------------------------------!
1816 SUBROUTINE chem_init_arrays
1817!
1818!-- Please use this place to allocate required arrays
1819
1820 END SUBROUTINE chem_init_arrays
1821
1822
1823!------------------------------------------------------------------------------!
1824! Description:
1825! ------------
1826!> Subroutine initializating chemistry_model_mod
1827!------------------------------------------------------------------------------!
1828 SUBROUTINE chem_init
1829
1830    USE chem_emissions_mod,                                                    &
1831        ONLY:  chem_emissions_init
1832       
1833    USE netcdf_data_input_mod,                                                 &
1834        ONLY:  init_3d
1835
1836
1837    INTEGER(iwp) ::  i !< running index x dimension
1838    INTEGER(iwp) ::  j !< running index y dimension
1839    INTEGER(iwp) ::  n !< running index for chemical species
1840
1841
1842    IF ( debug_output )  CALL debug_message( 'chem_init', 'start' )
1843!
1844!-- Next statement is to avoid compiler warning about unused variables
1845    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub + &
1846           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +    &
1847           ilu_urban ) == 0 )  CONTINUE
1848         
1849    IF ( emissions_anthropogenic )  CALL chem_emissions_init
1850!
1851!-- Chemistry variables will be initialized if availabe from dynamic
1852!-- input file. Note, it is possible to initialize only part of the chemistry
1853!-- variables from dynamic input.
1854    IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1855       DO  n = 1, nspec
1856          IF ( init_3d%from_file_chem(n) )  THEN
1857             DO  i = nxlg, nxrg
1858                DO  j = nysg, nyng
1859                   chem_species(n)%conc(:,j,i) = init_3d%chem_init(:,n)
1860                ENDDO
1861             ENDDO
1862          ENDIF
1863       ENDDO
1864    ENDIF
1865
1866    IF ( debug_output )  CALL debug_message( 'chem_init', 'end' )
1867
1868 END SUBROUTINE chem_init
1869
1870
1871!------------------------------------------------------------------------------!
1872! Description:
1873! ------------
1874!> Subroutine initializating chemistry_model_mod
1875!> internal workaround for chem_species dependency in chem_check_parameters
1876!------------------------------------------------------------------------------!
1877 SUBROUTINE chem_init_internal
1878 
1879    USE pegrid
1880
1881    USE netcdf_data_input_mod,                                                 &
1882        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,          &
1883               netcdf_data_input_chemistry_data
1884
1885!
1886!-- Local variables
1887    INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
1888    INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
1889    INTEGER(iwp) ::  lsp               !< running index for chem spcs
1890    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
1891
1892    IF ( emissions_anthropogenic )  THEN
1893       CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
1894    ENDIF
1895!
1896!-- Allocate memory for chemical species
1897    ALLOCATE( chem_species(nspec) )
1898    ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1899    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1900    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1901    ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) 
1902    ALLOCATE( phot_frequen(nphot) ) 
1903    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
1904    ALLOCATE( bc_cs_t_val(nspec) )
1905!
1906!-- Initialize arrays
1907    spec_conc_1 (:,:,:,:) = 0.0_wp
1908    spec_conc_2 (:,:,:,:) = 0.0_wp
1909    spec_conc_3 (:,:,:,:) = 0.0_wp
1910    spec_conc_av(:,:,:,:) = 0.0_wp
1911
1912
1913    DO  lsp = 1, nspec
1914       chem_species(lsp)%name    = spc_names(lsp)
1915
1916       chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
1917       chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
1918       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
1919       chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
1920
1921       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
1922       chem_species(lsp)%cssws_av    = 0.0_wp
1923!
1924!--    The following block can be useful when emission module is not applied. &
1925!--    if emission module is applied the following block will be overwritten.
1926       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1927       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1928       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) 
1929       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
1930       chem_species(lsp)%flux_s_cs = 0.0_wp                                     
1931       chem_species(lsp)%flux_l_cs = 0.0_wp                                     
1932       chem_species(lsp)%diss_s_cs = 0.0_wp                                     
1933       chem_species(lsp)%diss_l_cs = 0.0_wp                                     
1934!
1935!--   Allocate memory for initial concentration profiles
1936!--   (concentration values come from namelist)
1937!--   (@todo (FK): Because of this, chem_init is called in palm before
1938!--               check_parameters, since conc_pr_init is used there.
1939!--               We have to find another solution since chem_init should
1940!--               eventually be called from init_3d_model!!)
1941       ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
1942       chem_species(lsp)%conc_pr_init(:) = 0.0_wp
1943
1944    ENDDO
1945!
1946!-- Set control flags for decycling only at lateral boundary cores, within the
1947!-- inner cores the decycle flag is set to .False.. Even though it does not
1948!-- affect the setting of chemistry boundary conditions, this flag is used to
1949!-- set advection control flags appropriately.
1950    decycle_chem_lr = MERGE( decycle_chem_lr, .FALSE.,                         &
1951                             nxl == 0  .OR.  nxr == nx )
1952    decycle_chem_ns = MERGE( decycle_chem_ns, .FALSE.,                         &
1953                             nys == 0  .OR.  nyn == ny )
1954!
1955!-- For some passive scalars decycling may be enabled. This case, the lateral
1956!-- boundary conditions are non-cyclic for these scalars (chemical species
1957!-- and aerosols), while the other scalars may have
1958!-- cyclic boundary conditions. However, large gradients near the boundaries
1959!-- may produce stationary numerical oscillations near the lateral boundaries
1960!-- when a higher-order scheme is applied near these boundaries.
1961!-- To get rid-off this, set-up additional flags that control the order of the
1962!-- scalar advection scheme near the lateral boundaries for passive scalars
1963!-- with decycling.
1964    IF ( scalar_advec == 'ws-scheme' )  THEN
1965       ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1966!
1967!--    In case of decyling, set Neumann boundary conditions for wall_flags_0
1968!--    bit 31 instead of cyclic boundary conditions.
1969!--    Bit 31 is used to identify extended degradation zones (please see
1970!--    following comment).
1971!--    Note, since several also other modules like Salsa or other future
1972!--    one may access this bit but may have other boundary conditions, the
1973!--    original value of wall_flags_0 bit 31 must not be modified. Hence,
1974!--    store the boundary conditions directly on cs_advc_flags_s.
1975!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and
1976!--    bit 31 won't be used to control the numerical order.
1977!--    Initialize with flag 31 only.
1978       cs_advc_flags_s = 0
1979       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0,               &
1980                                BTEST( wall_flags_0, 31 ) )
1981
1982       IF ( decycle_chem_ns )  THEN
1983          IF ( nys == 0  )  THEN
1984             DO  i = 1, nbgp     
1985                cs_advc_flags_s(:,nys-i,:) = MERGE(                            &
1986                                        IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
1987                                        IBCLR( cs_advc_flags_s(:,nys,:), 31 ), &
1988                                        BTEST( cs_advc_flags_s(:,nys,:), 31 )  &
1989                                                  )
1990             ENDDO
1991          ENDIF
1992          IF ( nyn == ny )  THEN
1993             DO  i = 1, nbgp 
1994                cs_advc_flags_s(:,nyn+i,:) = MERGE(                            &
1995                                        IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
1996                                        IBCLR( cs_advc_flags_s(:,nyn,:), 31 ), &
1997                                        BTEST( cs_advc_flags_s(:,nyn,:), 31 )  &
1998                                                  )
1999             ENDDO
2000          ENDIF
2001       ENDIF
2002       IF ( .NOT. decycle_chem_lr )  THEN
2003          IF ( nxl == 0  )  THEN
2004             DO  i = 1, nbgp   
2005                cs_advc_flags_s(:,:,nxl-i) = MERGE(                            &
2006                                        IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
2007                                        IBCLR( cs_advc_flags_s(:,:,nxl), 31 ), &
2008                                        BTEST( cs_advc_flags_s(:,:,nxl), 31 )  &
2009                                                  )
2010             ENDDO
2011          ENDIF
2012          IF ( nxr == nx )  THEN
2013             DO  i = 1, nbgp   
2014                cs_advc_flags_s(:,:,nxr+i) = MERGE(                            &
2015                                        IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
2016                                        IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
2017                                        BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
2018                                                  )
2019             ENDDO
2020          ENDIF     
2021       ENDIF
2022!
2023!--    To initialize advection flags appropriately, pass the boundary flags.
2024!--    The last argument indicates that a passive scalar is treated, where
2025!--    the horizontal advection terms are degraded already 2 grid points before
2026!--    the lateral boundary to avoid stationary oscillations at large-gradients.
2027!--    Also, extended degradation zones are applied, where horizontal advection of
2028!--    passive scalars is discretized by first-order scheme at all grid points
2029!--    that in the vicinity of buildings (<= 3 grid points). Even though no
2030!--    building is within the numerical stencil, first-order scheme is used.
2031!--    At fourth and fifth grid point the order of the horizontal advection scheme
2032!--    is successively upgraded.
2033!--    These extended degradation zones are used to avoid stationary numerical
2034!--    oscillations, which are responsible for high concentration maxima that may
2035!--    appear under shear-free stable conditions.
2036       CALL ws_init_flags_scalar(                                              &
2037                  bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr, &
2038                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns, &
2039                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr, &
2040                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns, &
2041                  cs_advc_flags_s, .TRUE. )
2042    ENDIF
2043!
2044!-- Initial concentration of profiles is prescribed by parameters cs_profile
2045!-- and cs_heights in the namelist &chemistry_parameters
2046
2047    CALL chem_init_profiles
2048!   
2049!-- In case there is dynamic input file, create a list of names for chemistry
2050!-- initial input files. Also, initialize array that indicates whether the
2051!-- respective variable is on file or not.
2052    IF ( input_pids_dynamic )  THEN   
2053       ALLOCATE( init_3d%var_names_chem(1:nspec) )
2054       ALLOCATE( init_3d%from_file_chem(1:nspec) )
2055       init_3d%from_file_chem(:) = .FALSE.
2056       
2057       DO  lsp = 1, nspec
2058          init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
2059       ENDDO
2060    ENDIF
2061!
2062!-- Initialize model variables
2063    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2064         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2065!
2066!--    First model run of a possible job queue.
2067!--    Initial profiles of the variables must be computed.
2068       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2069!
2070!--       Transfer initial profiles to the arrays of the 3D model
2071          DO  lsp = 1, nspec
2072             DO  i = nxlg, nxrg
2073                DO  j = nysg, nyng
2074                   DO lpr_lev = 1, nz + 1
2075                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
2076                   ENDDO
2077                ENDDO
2078             ENDDO   
2079          ENDDO
2080
2081       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
2082       THEN
2083
2084          DO  lsp = 1, nspec 
2085             DO  i = nxlg, nxrg
2086                DO  j = nysg, nyng
2087                   chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init   
2088                ENDDO
2089             ENDDO
2090          ENDDO
2091
2092       ENDIF
2093!
2094!--    If required, change the surface chem spcs at the start of the 3D run
2095       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN           
2096          DO  lsp = 1, nspec 
2097             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
2098                                               cs_surface_initial_change(lsp)
2099          ENDDO
2100       ENDIF 
2101!
2102!--    Initiale old and new time levels.
2103       DO  lsp = 1, nvar
2104          chem_species(lsp)%tconc_m = 0.0_wp                     
2105          chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
2106       ENDDO
2107
2108    ENDIF
2109
2110    DO  lsp = 1, nphot
2111       phot_frequen(lsp)%name = phot_names(lsp)
2112!
2113!-- todo: remove or replace by "CALL message" mechanism (kanani)
2114!--       IF( myid == 0 )  THEN
2115!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
2116!--       ENDIF
2117          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
2118    ENDDO
2119
2120!    CALL photolysis_init   ! probably also required for restart
2121
2122    RETURN
2123
2124 END SUBROUTINE chem_init_internal
2125
2126
2127!------------------------------------------------------------------------------!
2128! Description:
2129! ------------
2130!> Subroutine defining initial vertical profiles of chemical species (given by
2131!> namelist parameters chem_profiles and chem_heights)  --> which should work
2132!> analogue to parameters u_profile, v_profile and uv_heights)
2133!------------------------------------------------------------------------------!
2134 SUBROUTINE chem_init_profiles             
2135!
2136!-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
2137!< We still need to see what has to be done in case of restart run
2138
2139    USE chem_modules
2140
2141!
2142!-- Local variables
2143    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
2144    INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
2145    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
2146    INTEGER ::  npr_lev    !< the next available profile lev
2147!
2148!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
2149!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
2150!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
2151!-- "cs_heights" are prescribed, their values will!override the constant profile given by
2152!-- "cs_surface".
2153    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2154       lsp_usr = 1
2155       DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
2156          DO  lsp = 1, nspec                                !
2157!
2158!--          create initial profile (conc_pr_init) for each chemical species
2159             IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
2160                IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
2161!
2162!--               set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
2163                   DO lpr_lev = 0, nzt+1
2164                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
2165                   ENDDO
2166                ELSE
2167                   IF ( cs_heights(1,1) /= 0.0_wp )  THEN
2168                      message_string = 'The surface value of cs_heights must be 0.0'
2169                      CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
2170                   ENDIF
2171
2172                   use_prescribed_profile_data = .TRUE.
2173
2174                   npr_lev = 1
2175!                   chem_species(lsp)%conc_pr_init(0) = 0.0_wp
2176                   DO  lpr_lev = 1, nz+1
2177                      IF ( npr_lev < 100 )  THEN
2178                         DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
2179                            npr_lev = npr_lev + 1
2180                            IF ( npr_lev == 100 )  THEN
2181                               message_string = 'number of chem spcs exceeding the limit'
2182                               CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
2183                               EXIT
2184                            ENDIF
2185                         ENDDO
2186                      ENDIF
2187                      IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
2188                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
2189                              ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
2190                              ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
2191                              ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
2192                      ELSE
2193                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
2194                      ENDIF
2195                   ENDDO
2196                ENDIF
2197!
2198!--          If a profile is prescribed explicity using cs_profiles and cs_heights, then 
2199!--          chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
2200!--          on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
2201             ENDIF
2202          ENDDO
2203          lsp_usr = lsp_usr + 1
2204       ENDDO
2205    ENDIF
2206
2207 END SUBROUTINE chem_init_profiles
2208
2209 
2210!------------------------------------------------------------------------------!
2211! Description:
2212! ------------
2213!> Subroutine to integrate chemical species in the given chemical mechanism
2214!------------------------------------------------------------------------------!
2215 SUBROUTINE chem_integrate_ij( i, j )
2216
2217    USE statistics,                                                          &
2218        ONLY:  weight_pres
2219
2220    USE control_parameters,                                                  &
2221        ONLY:  dt_3d, intermediate_timestep_count, time_since_reference_point
2222
2223
2224    INTEGER,INTENT(IN)       :: i
2225    INTEGER,INTENT(IN)       :: j
2226!
2227!--   local variables
2228    INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
2229    INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
2230    INTEGER, DIMENSION(20)    :: istatus
2231    REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
2232    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
2233    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
2234    REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
2235    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
2236    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
2237                                                                              !<    molecules cm^{-3} and ppm
2238
2239    INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
2240    INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
2241
2242    REAL(wp)                         ::  conv                                !< conversion factor
2243    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
2244    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp                  !< Conversion factor fraction to ppm
2245!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
2246!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
2247    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
2248    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
2249    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^3
2250    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
2251
2252    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
2253
2254    REAL(kind=wp)  :: dt_chem                                             
2255!
2256!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
2257    IF (chem_gasphase_on)  THEN
2258       nacc = 0
2259       nrej = 0
2260
2261       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
2262!
2263!--    convert ppm to molecules/cm**3
2264!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
2265!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
2266       conv = ppm2fr * xna / vmolcm
2267       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
2268       tmp_fact_i = 1.0_wp/tmp_fact
2269
2270       IF ( humidity )  THEN
2271          IF ( bulk_cloud_model )  THEN
2272             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
2273                             xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2274          ELSE
2275             tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2276          ENDIF
2277       ELSE
2278          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)          !< Constant value for q if water vapor is not computed
2279       ENDIF
2280
2281       DO  lsp = 1,nspec
2282          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
2283       ENDDO
2284
2285       DO lph = 1,nphot
2286          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
2287       ENDDO
2288!
2289!--    Compute length of time step
2290       IF ( call_chem_at_all_substeps )  THEN
2291          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
2292       ELSE
2293          dt_chem = dt_3d
2294       ENDIF
2295
2296       cs_time_step = dt_chem
2297
2298       IF(maxval(rcntrl) > 0.0)  THEN    ! Only if rcntrl is set
2299          IF( time_since_reference_point <= 2*dt_3d)  THEN
2300             rcntrl_local = 0
2301          ELSE
2302             rcntrl_local = rcntrl
2303          ENDIF
2304       ELSE
2305          rcntrl_local = 0
2306       END IF
2307
2308       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
2309            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
2310
2311       DO  lsp = 1,nspec
2312          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
2313       ENDDO
2314
2315
2316    ENDIF
2317
2318    RETURN
2319 END SUBROUTINE chem_integrate_ij
2320
2321
2322!------------------------------------------------------------------------------!
2323! Description:
2324! ------------
2325!> Subroutine defining parin for &chemistry_parameters for chemistry model
2326!------------------------------------------------------------------------------!
2327 SUBROUTINE chem_parin
2328
2329    USE chem_modules
2330    USE control_parameters
2331
2332    USE pegrid
2333    USE statistics
2334
2335
2336    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
2337
2338    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
2339    INTEGER(iwp) ::  i                                 !<
2340    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
2341
2342
2343    NAMELIST /chemistry_parameters/  bc_cs_b,                          &
2344         bc_cs_t,                          &
2345         call_chem_at_all_substeps,        &
2346         chem_debug0,                      &
2347         chem_debug1,                      &
2348         chem_debug2,                      &
2349         chem_gasphase_on,                 &
2350         chem_mechanism,                   &         
2351         cs_heights,                       &
2352         cs_name,                          &
2353         cs_profile,                       &
2354         cs_surface,                       &
2355         cs_surface_initial_change,        &
2356         cs_vertical_gradient_level,       &
2357         daytype_mdh,                      &
2358         decycle_chem_lr,                  &
2359         decycle_chem_ns,                  &           
2360         decycle_method,                   &
2361         deposition_dry,                   &
2362         emissions_anthropogenic,          &
2363         emiss_lod,                        &
2364         emiss_factor_main,                &
2365         emiss_factor_side,                &                     
2366         icntrl,                           &
2367         main_street_id,                   &
2368         max_street_id,                    &
2369         mode_emis,                        &
2370         my_steps,                         &
2371         rcntrl,                           &
2372         side_street_id,                   &
2373         photolysis_scheme,                &
2374         wall_csflux,                      &
2375         cs_vertical_gradient,             &
2376         top_csflux,                       & 
2377         surface_csflux,                   &
2378         surface_csflux_name,              &
2379         time_fac_type
2380!
2381!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
2382!-- so this way we could prescribe a specific flux value for each species
2383    !>  chemistry_parameters for initial profiles
2384    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
2385    !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
2386    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
2387    !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
2388    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
2389    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
2390!
2391!--   Read chem namelist   
2392
2393    CHARACTER(LEN=8)    :: solver_type
2394
2395    icntrl    = 0
2396    rcntrl    = 0.0_wp
2397    my_steps  = 0.0_wp
2398    photolysis_scheme = 'simple'
2399    atol = 1.0_wp
2400    rtol = 0.01_wp
2401!
2402!--   Try to find chemistry package
2403    REWIND ( 11 )
2404    line = ' '
2405    DO   WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
2406       READ ( 11, '(A)', END=20 )  line
2407    ENDDO
2408    BACKSPACE ( 11 )
2409!
2410!--   Read chemistry namelist
2411    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
2412!
2413!--   Enable chemistry model
2414    air_chemistry = .TRUE.                   
2415    GOTO 20
2416
2417 10 BACKSPACE( 11 )
2418    READ( 11 , '(A)') line
2419    CALL parin_fail_message( 'chemistry_parameters', line )
2420
2421 20 CONTINUE
2422
2423!
2424!-- synchronize emiss_lod and mod_emis only if emissions_anthropogenic
2425!-- is activated in the namelist.  Otherwise their values are "don't care"
2426    IF ( emissions_anthropogenic )  THEN
2427!
2428!--    check for emission mode for chem species
2429
2430       IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
2431          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.      &
2432               ( mode_emis /= 'DEFAULT'        )    .AND.      &
2433               ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2434             message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2435             CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2436          ENDIF
2437       ELSE
2438          IF ( ( emiss_lod /= 0 )    .AND.         &
2439               ( emiss_lod /= 1 )    .AND.         &
2440               ( emiss_lod /= 2 ) )  THEN
2441             message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
2442             CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2443          ENDIF
2444       ENDIF
2445
2446!
2447! for reference (ecc)
2448!    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2449!       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2450!       CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2451!    ENDIF
2452
2453!
2454!-- conflict resolution for emiss_lod and mode_emis
2455!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
2456!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
2457!-- this check is in place to retain backward compatibility with salsa until the
2458!-- code is migrated completed to emiss_lod
2459!-- note that
2460
2461       IF  ( emiss_lod >= 0 ) THEN
2462
2463          SELECT CASE  ( emiss_lod )
2464             CASE (0)  !- parameterized mode
2465                mode_emis = 'PARAMETERIZED'
2466             CASE (1)  !- default mode
2467                mode_emis = 'DEFAULT'
2468             CASE (2)  !- preprocessed mode
2469                mode_emis = 'PRE-PROCESSED'
2470          END SELECT
2471       
2472          message_string = 'Synchronizing mode_emis to defined emiss_lod'               //  &
2473                           CHAR(10)  //  '          '                                   //  &
2474                           'NOTE - mode_emis will be depreciated in future releases'    //  &
2475                           CHAR(10)  //  '          '                                   //  &
2476                           'please use emiss_lod to define emission mode'
2477 
2478          CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
2479
2480       ELSE ! if emiss_lod is not set
2481
2482          SELECT CASE ( mode_emis )
2483             CASE ('PARAMETERIZED')
2484                emiss_lod = 0
2485             CASE ('DEFAULT')
2486                emiss_lod = 1
2487             CASE ('PRE-PROCESSED')
2488                emiss_lod = 2
2489          END SELECT
2490
2491          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //  &
2492                           CHAR(10)  //  '          '                                   //  &
2493                           'NOTE - mode_emis will be depreciated in future releases'    //  &
2494                           CHAR(10)  //  '          '                                   //  &
2495                           '       please use emiss_lod to define emission mode'
2496
2497          CALL message ( 'parin_chem', 'CM0464', 0, 0, 0, 6, 0 )
2498       ENDIF
2499
2500    ENDIF  ! if emissions_anthropengic
2501
2502    t_steps = my_steps         
2503!
2504!--    Determine the number of user-defined profiles and append them to the
2505!--    standard data output (data_output_pr)
2506    max_pr_cs_tmp = 0 
2507    i = 1
2508
2509    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
2510       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
2511          max_pr_cs_tmp = max_pr_cs_tmp+1
2512       ENDIF
2513       i = i +1
2514    ENDDO
2515
2516    IF ( max_pr_cs_tmp > 0 )  THEN
2517       cs_pr_namelist_found = .TRUE.
2518       max_pr_cs = max_pr_cs_tmp
2519    ENDIF
2520
2521    !      Set Solver Type
2522    IF(icntrl(3) == 0)  THEN
2523       solver_type = 'rodas3'           !Default
2524    ELSE IF(icntrl(3) == 1)  THEN
2525       solver_type = 'ros2'
2526    ELSE IF(icntrl(3) == 2)  THEN
2527       solver_type = 'ros3'
2528    ELSE IF(icntrl(3) == 3)  THEN
2529       solver_type = 'ro4'
2530    ELSE IF(icntrl(3) == 4)  THEN
2531       solver_type = 'rodas3'
2532    ELSE IF(icntrl(3) == 5)  THEN
2533       solver_type = 'rodas4'
2534    ELSE IF(icntrl(3) == 6)  THEN
2535       solver_type = 'Rang3'
2536    ELSE
2537       message_string = 'illegal solver type'
2538       CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
2539    END IF
2540
2541!
2542!--   todo: remove or replace by "CALL message" mechanism (kanani)
2543!       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
2544!kk    Has to be changed to right calling sequence
2545!        IF(myid == 0)  THEN
2546!           write(9,*) ' '
2547!           write(9,*) 'kpp setup '
2548!           write(9,*) ' '
2549!           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM( solver_type )
2550!           write(9,*) ' '
2551!           write(9,*) '    Hstart  = ',rcntrl(3)
2552!           write(9,*) '    FacMin  = ',rcntrl(4)
2553!           write(9,*) '    FacMax  = ',rcntrl(5)
2554!           write(9,*) ' '
2555!           IF(vl_dim > 1)  THEN
2556!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
2557!           ELSE
2558!              write(9,*) '    Scalar mode'
2559!           ENDIF
2560!           write(9,*) ' '
2561!        END IF
2562
2563    RETURN
2564
2565 END SUBROUTINE chem_parin
2566
2567
2568!------------------------------------------------------------------------------!
2569! Description:
2570! ------------
2571!> Call for all grid points
2572!------------------------------------------------------------------------------!
2573    SUBROUTINE chem_actions( location )
2574
2575
2576    CHARACTER (LEN=*), INTENT(IN) ::  location !< call location string
2577
2578    SELECT CASE ( location )
2579
2580       CASE ( 'before_prognostic_equations' )
2581!
2582!--       Chemical reactions and deposition
2583          IF ( chem_gasphase_on ) THEN
2584!
2585!--          If required, calculate photolysis frequencies -
2586!--          UNFINISHED: Why not before the intermediate timestep loop?
2587             IF ( intermediate_timestep_count ==  1 )  THEN
2588                CALL photolysis_control
2589             ENDIF
2590
2591          ENDIF
2592
2593       CASE DEFAULT
2594          CONTINUE
2595
2596    END SELECT
2597
2598    END SUBROUTINE chem_actions
2599
2600
2601!------------------------------------------------------------------------------!
2602! Description:
2603! ------------
2604!> Call for grid points i,j
2605!------------------------------------------------------------------------------!
2606
2607    SUBROUTINE chem_actions_ij( i, j, location )
2608
2609
2610    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
2611    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
2612    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
2613    INTEGER(iwp)  ::  dummy  !< call location string
2614
2615    IF ( air_chemistry    )   dummy = i + j
2616
2617    SELECT CASE ( location )
2618
2619       CASE DEFAULT
2620          CONTINUE
2621
2622    END SELECT
2623
2624
2625    END SUBROUTINE chem_actions_ij
2626
2627
2628!------------------------------------------------------------------------------!
2629! Description:
2630! ------------
2631!> Call for all grid points
2632!------------------------------------------------------------------------------!
2633    SUBROUTINE chem_non_advective_processes()
2634
2635
2636      INTEGER(iwp) ::  i  !<
2637      INTEGER(iwp) ::  j  !<
2638
2639      !
2640      !-- Calculation of chemical reactions and deposition.
2641
2642
2643      IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
2644
2645         IF ( chem_gasphase_on ) THEN
2646            CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2647            !$OMP PARALLEL PRIVATE (i,j)
2648            !$OMP DO schedule(static,1)
2649            DO  i = nxl, nxr
2650               DO  j = nys, nyn
2651                  CALL chem_integrate( i, j )
2652               ENDDO
2653            ENDDO
2654            !$OMP END PARALLEL
2655            CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2656         ENDIF
2657
2658         IF ( deposition_dry )  THEN
2659            CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2660            DO  i = nxl, nxr
2661               DO  j = nys, nyn
2662                  CALL chem_depo( i, j )
2663               ENDDO
2664            ENDDO
2665            CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2666         ENDIF
2667
2668      ENDIF
2669
2670
2671
2672    END SUBROUTINE chem_non_advective_processes
2673
2674
2675!------------------------------------------------------------------------------!
2676! Description:
2677! ------------
2678!> Call for grid points i,j
2679!------------------------------------------------------------------------------!
2680 SUBROUTINE chem_non_advective_processes_ij( i, j )
2681
2682
2683   INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
2684   INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
2685
2686!
2687!-- Calculation of chemical reactions and deposition.
2688
2689
2690   IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
2691
2692      IF ( chem_gasphase_on ) THEN
2693         CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2694         CALL chem_integrate( i, j )
2695         CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2696      ENDIF
2697
2698      IF ( deposition_dry )  THEN
2699         CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2700         CALL chem_depo( i, j )
2701         CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2702      ENDIF
2703
2704   ENDIF
2705
2706
2707
2708 END SUBROUTINE chem_non_advective_processes_ij
2709 
2710!------------------------------------------------------------------------------!
2711! Description:
2712! ------------
2713!> routine for exchange horiz of chemical quantities 
2714!------------------------------------------------------------------------------!
2715 SUBROUTINE chem_exchange_horiz_bounds
2716 
2717   INTEGER(iwp) ::  lsp       !<
2718   INTEGER(iwp) ::  lsp_usr   !<
2719 
2720!
2721!--    Loop over chemical species       
2722       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
2723       DO  lsp = 1, nvar
2724          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
2725          lsp_usr = 1 
2726          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
2727             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
2728
2729                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                     &
2730                                                  chem_species(lsp)%conc_pr_init ) 
2731             
2732             ENDIF
2733             lsp_usr = lsp_usr +1
2734          ENDDO
2735
2736
2737       ENDDO
2738       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
2739
2740
2741 END SUBROUTINE chem_exchange_horiz_bounds
2742
2743 
2744!------------------------------------------------------------------------------!
2745! Description:
2746! ------------
2747!> Subroutine calculating prognostic equations for chemical species
2748!> (vector-optimized).
2749!> Routine is called separately for each chemical species over a loop from
2750!> prognostic_equations.
2751!------------------------------------------------------------------------------!
2752 SUBROUTINE chem_prognostic_equations()
2753
2754
2755    INTEGER ::  i   !< running index
2756    INTEGER ::  j   !< running index
2757    INTEGER ::  k   !< running index
2758
2759    INTEGER(iwp) ::  ilsp   !<
2760
2761
2762    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2763
2764    DO  ilsp = 1, nvar
2765!
2766!--    Tendency terms for chemical species
2767       tend = 0.0_wp
2768!
2769!--    Advection terms
2770       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2771          IF ( ws_scheme_sca )  THEN
2772             CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
2773                              bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr,         &
2774                              bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns,         &
2775                              bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr,         &
2776                              bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns )
2777          ELSE
2778             CALL advec_s_pw( chem_species(ilsp)%conc )
2779          ENDIF
2780       ELSE
2781          CALL advec_s_up( chem_species(ilsp)%conc )
2782       ENDIF
2783!
2784!--    Diffusion terms  (the last three arguments are zero)
2785       CALL diffusion_s( chem_species(ilsp)%conc,                                                  &
2786            surf_def_h(0)%cssws(ilsp,:),                                                           &
2787            surf_def_h(1)%cssws(ilsp,:),                                                           &
2788            surf_def_h(2)%cssws(ilsp,:),                                                           &
2789            surf_lsm_h%cssws(ilsp,:),                                                              &
2790            surf_usm_h%cssws(ilsp,:),                                                              &
2791            surf_def_v(0)%cssws(ilsp,:),                                                           &
2792            surf_def_v(1)%cssws(ilsp,:),                                                           &
2793            surf_def_v(2)%cssws(ilsp,:),                                                           &
2794            surf_def_v(3)%cssws(ilsp,:),                                                           &
2795            surf_lsm_v(0)%cssws(ilsp,:),                                                           &
2796            surf_lsm_v(1)%cssws(ilsp,:),                                                           &
2797            surf_lsm_v(2)%cssws(ilsp,:),                                                           &
2798            surf_lsm_v(3)%cssws(ilsp,:),                                                           &
2799            surf_usm_v(0)%cssws(ilsp,:),                                                           &
2800            surf_usm_v(1)%cssws(ilsp,:),                                                           &
2801            surf_usm_v(2)%cssws(ilsp,:),                                                           &
2802            surf_usm_v(3)%cssws(ilsp,:) )
2803!
2804!--    Prognostic equation for chemical species
2805       DO  i = nxl, nxr
2806          DO  j = nys, nyn
2807             DO  k = nzb+1, nzt
2808                chem_species(ilsp)%conc_p(k,j,i) =   chem_species(ilsp)%conc(k,j,i)                &
2809                     + ( dt_3d  *                                                                  &
2810                     (   tsc(2) * tend(k,j,i)                                                      &
2811                     + tsc(3) * chem_species(ilsp)%tconc_m(k,j,i)                                  &
2812                     )                                                                             &
2813                     - tsc(5) * rdf_sc(k)                                                          &
2814                     * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )     &
2815                     )                                                                             &
2816                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2817
2818                IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
2819                   chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)
2820                ENDIF
2821             ENDDO
2822          ENDDO
2823       ENDDO
2824!
2825!--    Calculate tendencies for the next Runge-Kutta step
2826       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2827          IF ( intermediate_timestep_count == 1 )  THEN
2828             DO  i = nxl, nxr
2829                DO  j = nys, nyn
2830                   DO  k = nzb+1, nzt
2831                      chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2832                   ENDDO
2833                ENDDO
2834             ENDDO
2835          ELSEIF ( intermediate_timestep_count < &
2836               intermediate_timestep_count_max )  THEN
2837             DO  i = nxl, nxr
2838                DO  j = nys, nyn
2839                   DO  k = nzb+1, nzt
2840                      chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
2841                           + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2842                   ENDDO
2843                ENDDO
2844             ENDDO
2845          ENDIF
2846       ENDIF
2847
2848    ENDDO
2849
2850    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
2851
2852 END SUBROUTINE chem_prognostic_equations
2853
2854
2855!------------------------------------------------------------------------------!
2856! Description:
2857! ------------
2858!> Subroutine calculating prognostic equations for chemical species
2859!> (cache-optimized).
2860!> Routine is called separately for each chemical species over a loop from
2861!> prognostic_equations.
2862!------------------------------------------------------------------------------!
2863 SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
2864
2865
2866    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
2867    INTEGER(iwp) :: ilsp
2868!
2869!-- local variables
2870
2871    INTEGER :: k
2872
2873    DO  ilsp = 1, nvar
2874!
2875!--    Tendency-terms for chem spcs.
2876       tend(:,j,i) = 0.0_wp
2877!
2878!--    Advection terms
2879       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2880          IF ( ws_scheme_sca )  THEN
2881             CALL advec_s_ws( cs_advc_flags_s,                                                     &
2882                              i,                                                                   &
2883                              j,                                                                   &
2884                              chem_species(ilsp)%conc,                                             &
2885                              'kc',                                                                &
2886                              chem_species(ilsp)%flux_s_cs,                                        &
2887                              chem_species(ilsp)%diss_s_cs,                                        &
2888                              chem_species(ilsp)%flux_l_cs,                                        &
2889                              chem_species(ilsp)%diss_l_cs,                                        &
2890                              i_omp_start,                                                         &
2891                              tn,                                                                  &
2892                              bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr,         &
2893                              bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns,         &
2894                              bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr,         &
2895                              bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns,         &
2896                              monotonic_limiter_z )
2897          ELSE
2898             CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
2899          ENDIF
2900       ELSE
2901          CALL advec_s_up( i, j, chem_species(ilsp)%conc )
2902       ENDIF
2903!
2904!--    Diffusion terms (the last three arguments are zero)
2905
2906       CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
2907            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
2908            surf_def_h(2)%cssws(ilsp,:),                                                           &
2909            surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),                                    &
2910            surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),                              &
2911            surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),                              &
2912            surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),                              &
2913            surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),                              &
2914            surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),                              &
2915            surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2916!
2917!--    Prognostic equation for chem spcs
2918       DO  k = nzb+1, nzt
2919          chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) + ( dt_3d  *           &
2920               ( tsc(2) * tend(k,j,i) +                                                            &
2921               tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) )                                        &
2922               - tsc(5) * rdf_sc(k)                                                                &
2923               * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
2924               )                                                                                   &
2925               * MERGE( 1.0_wp, 0.0_wp,                                                            &
2926               BTEST( wall_flags_0(k,j,i), 0 )                                                     &
2927               )
2928
2929          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN               
2930             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6               
2931          ENDIF
2932       ENDDO
2933!
2934!--    Calculate tendencies for the next Runge-Kutta step
2935       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2936          IF ( intermediate_timestep_count == 1 )  THEN
2937             DO  k = nzb+1, nzt
2938                chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2939             ENDDO
2940          ELSEIF ( intermediate_timestep_count < &
2941               intermediate_timestep_count_max )  THEN
2942             DO  k = nzb+1, nzt
2943                chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
2944                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2945             ENDDO
2946          ENDIF
2947       ENDIF
2948
2949    ENDDO
2950
2951 END SUBROUTINE chem_prognostic_equations_ij
2952
2953
2954!------------------------------------------------------------------------------!
2955! Description:
2956! ------------
2957!> Subroutine to read restart data of chemical species
2958!------------------------------------------------------------------------------!
2959 SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
2960                            nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
2961                            nys_on_file, tmp_3d, found )
2962
2963    USE control_parameters
2964
2965
2966    CHARACTER (LEN=20) :: spc_name_av !<   
2967
2968    INTEGER(iwp) ::  lsp             !<
2969    INTEGER(iwp) ::  k               !<
2970    INTEGER(iwp) ::  nxlc            !<
2971    INTEGER(iwp) ::  nxlf            !<
2972    INTEGER(iwp) ::  nxl_on_file     !<   
2973    INTEGER(iwp) ::  nxrc            !<
2974    INTEGER(iwp) ::  nxrf            !<
2975    INTEGER(iwp) ::  nxr_on_file     !<   
2976    INTEGER(iwp) ::  nync            !<
2977    INTEGER(iwp) ::  nynf            !<
2978    INTEGER(iwp) ::  nyn_on_file     !<   
2979    INTEGER(iwp) ::  nysc            !<
2980    INTEGER(iwp) ::  nysf            !<
2981    INTEGER(iwp) ::  nys_on_file     !<   
2982
2983    LOGICAL, INTENT(OUT) :: found 
2984
2985    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
2986
2987
2988    found = .FALSE. 
2989
2990
2991    IF ( ALLOCATED(chem_species) )  THEN
2992
2993       DO  lsp = 1, nspec
2994
2995          !< for time-averaged chemical conc.
2996          spc_name_av  =  TRIM( chem_species(lsp)%name )//'_av'
2997
2998          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )    &
2999             THEN
3000             !< read data into tmp_3d
3001             IF ( k == 1 )  READ ( 13 )  tmp_3d 
3002             !< fill ..%conc in the restart run   
3003             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                    &
3004                  nxlc-nbgp:nxrc+nbgp) =                  & 
3005                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3006             found = .TRUE.
3007          ELSEIF (restart_string(1:length) == spc_name_av )  THEN
3008             IF ( k == 1 )  READ ( 13 )  tmp_3d
3009             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,                 &
3010                  nxlc-nbgp:nxrc+nbgp) =               &
3011                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3012             found = .TRUE.
3013          ENDIF
3014
3015       ENDDO
3016
3017    ENDIF
3018
3019
3020 END SUBROUTINE chem_rrd_local
3021
3022
3023!-------------------------------------------------------------------------------!
3024!> Description:
3025!> Calculation of horizontally averaged profiles
3026!> This routine is called for every statistic region (sr) defined by the user,
3027!> but at least for the region "total domain" (sr=0).
3028!> quantities.
3029!-------------------------------------------------------------------------------!
3030 SUBROUTINE chem_statistics( mode, sr, tn )
3031
3032
3033    USE arrays_3d
3034
3035    USE statistics
3036
3037
3038    CHARACTER (LEN=*) ::  mode   !<
3039
3040    INTEGER(iwp) ::  i    !< running index on x-axis
3041    INTEGER(iwp) ::  j    !< running index on y-axis
3042    INTEGER(iwp) ::  k    !< vertical index counter
3043    INTEGER(iwp) ::  sr   !< statistical region
3044    INTEGER(iwp) ::  tn   !< thread number
3045    INTEGER(iwp) ::  lpr  !< running index chem spcs
3046
3047    IF ( mode == 'profiles' )  THEN
3048       !
3049!
3050!--    Sample on how to calculate horizontally averaged profiles of user-
3051!--    defined quantities. Each quantity is identified by the index
3052!--    "pr_palm+#" where "#" is an integer starting from 1. These
3053!--    user-profile-numbers must also be assigned to the respective strings
3054!--    given by data_output_pr_cs in routine user_check_data_output_pr.
3055!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
3056!--                     w*pt*), dim-4 = statistical region.
3057
3058!$OMP DO
3059       DO  i = nxl, nxr
3060          DO  j = nys, nyn
3061             DO  k = nzb, nzt+1
3062                DO lpr = 1, cs_pr_count
3063
3064                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
3065                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
3066                        rmask(j,i,sr)  *                                   &
3067                        MERGE( 1.0_wp, 0.0_wp,                             &
3068                        BTEST( wall_flags_0(k,j,i), 22 ) )
3069                ENDDO
3070             ENDDO
3071          ENDDO
3072       ENDDO
3073    ELSEIF ( mode == 'time_series' )  THEN
3074!      @todo
3075    ENDIF
3076
3077 END SUBROUTINE chem_statistics
3078
3079
3080!------------------------------------------------------------------------------!
3081! Description:
3082! ------------
3083!> Subroutine for swapping of timelevels for chemical species
3084!> called out from subroutine swap_timelevel
3085!------------------------------------------------------------------------------!
3086
3087
3088 SUBROUTINE chem_swap_timelevel( level )
3089
3090
3091    INTEGER(iwp), INTENT(IN) ::  level
3092!
3093!-- local variables
3094    INTEGER(iwp)             ::  lsp
3095
3096
3097    IF ( level == 0 )  THEN
3098       DO  lsp=1, nvar                                       
3099          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
3100          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
3101       ENDDO
3102    ELSE
3103       DO  lsp=1, nvar                                       
3104          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
3105          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
3106       ENDDO
3107    ENDIF
3108
3109    RETURN
3110 END SUBROUTINE chem_swap_timelevel
3111
3112
3113!------------------------------------------------------------------------------!
3114! Description:
3115! ------------
3116!> Subroutine to write restart data for chemistry model
3117!------------------------------------------------------------------------------!
3118 SUBROUTINE chem_wrd_local
3119
3120
3121    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
3122
3123    DO  lsp = 1, nspec
3124       CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
3125       WRITE ( 14 )  chem_species(lsp)%conc
3126       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
3127       WRITE ( 14 )  chem_species(lsp)%conc_av
3128    ENDDO
3129
3130 END SUBROUTINE chem_wrd_local
3131
3132
3133!-------------------------------------------------------------------------------!
3134! Description:
3135! ------------
3136!> Subroutine to calculate the deposition of gases and PMs. For now deposition
3137!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
3138!> considered. The deposition of particles is derived following Zhang et al.,
3139!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
3140!>     
3141!> @TODO: Consider deposition on vertical surfaces   
3142!> @TODO: Consider overlaying horizontal surfaces
3143!> @TODO: Consider resolved vegetation   
3144!> @TODO: Check error messages
3145!-------------------------------------------------------------------------------!
3146 SUBROUTINE chem_depo( i, j )
3147
3148    USE control_parameters,                                                 &   
3149         ONLY:  dt_3d, intermediate_timestep_count, latitude
3150
3151    USE arrays_3d,                                                          &
3152         ONLY:  dzw, rho_air_zw
3153
3154    USE date_and_time_mod,                                                  &
3155         ONLY:  day_of_year
3156
3157    USE surface_mod,                                                        &
3158         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
3159         surf_type, surf_usm_h
3160
3161    USE radiation_model_mod,                                                &
3162         ONLY:  cos_zenith
3163
3164
3165    INTEGER(iwp), INTENT(IN) ::  i
3166    INTEGER(iwp), INTENT(IN) ::  j
3167    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
3168    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
3169    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current surface element
3170    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
3171    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
3172    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
3173    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
3174    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
3175    INTEGER(iwp) ::  luv_dep                       !< matching DEPAC LU to luv_palm
3176    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
3177    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
3178    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
3179    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
3180    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
3181    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
3182
3183    INTEGER(iwp) ::  pspec                         !< running index
3184    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
3185!
3186!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
3187    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
3188    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
3189    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
3190    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
3191    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
3192    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
3193    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
3194    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
3195    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
3196    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
3197    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
3198    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
3199    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
3200    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
3201    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
3202    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
3203    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
3204    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
3205    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
3206!
3207!-- Water
3208    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM 
3209    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
3210    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
3211    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
3212    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
3213    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
3214!
3215!-- Pavement
3216    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
3217    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
3218    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
3219    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
3220    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
3221    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
3222    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
3223    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
3224    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
3225    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
3226    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
3227    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
3228    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
3229    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
3230    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
3231    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
3232!
3233!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
3234    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
3235    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
3236    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
3237
3238    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
3239
3240    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
3241
3242    REAL(wp) ::  dt_chem                !< length of chem time step
3243    REAL(wp) ::  dh                     !< vertical grid size
3244    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size
3245    REAL(wp) ::  dt_dh                  !< dt_chem/dh
3246
3247    REAL(wp) ::  dens              !< density at layer k at i,j 
3248    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
3249    REAL(wp) ::  ustar_surf        !< ustar at current surface element
3250    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
3251    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element
3252    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
3253    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
3254    REAL(wp) ::  lai               !< leaf area index at current surface element
3255    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
3256
3257    REAL(wp) ::  slinnfac       
3258    REAL(wp) ::  visc              !< Viscosity
3259    REAL(wp) ::  vs                !< Sedimentation velocity
3260    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
3261    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
3262    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
3263    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
3264
3265    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
3266    REAL(wp) ::  diffusivity       !< diffusivity
3267
3268
3269    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
3270    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
3271    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
3272    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
3273    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
3274    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
3275    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
3276    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
3277    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
3278                                                 
3279
3280    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
3281    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
3282    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
3283!
3284!-- Particle parameters (PM10 (1), PM25 (2))
3285!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
3286!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
3287    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
3288         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
3289         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
3290         /), (/ 3, 2 /) )
3291
3292    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
3293    LOGICAL ::  match_usm     !< flag indicating urban-type surface
3294!
3295!-- List of names of possible tracers
3296    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
3297         'NO2           ', &    !< NO2
3298         'NO            ', &    !< NO
3299         'O3            ', &    !< O3
3300         'CO            ', &    !< CO
3301         'form          ', &    !< FORM
3302         'ald           ', &    !< ALD
3303         'pan           ', &    !< PAN
3304         'mgly          ', &    !< MGLY
3305         'par           ', &    !< PAR
3306         'ole           ', &    !< OLE
3307         'eth           ', &    !< ETH
3308         'tol           ', &    !< TOL
3309         'cres          ', &    !< CRES
3310         'xyl           ', &    !< XYL
3311         'SO4a_f        ', &    !< SO4a_f
3312         'SO2           ', &    !< SO2
3313         'HNO2          ', &    !< HNO2
3314         'CH4           ', &    !< CH4
3315         'NH3           ', &    !< NH3
3316         'NO3           ', &    !< NO3
3317         'OH            ', &    !< OH
3318         'HO2           ', &    !< HO2
3319         'N2O5          ', &    !< N2O5
3320         'SO4a_c        ', &    !< SO4a_c
3321         'NH4a_f        ', &    !< NH4a_f
3322         'NO3a_f        ', &    !< NO3a_f
3323         'NO3a_c        ', &    !< NO3a_c
3324         'C2O3          ', &    !< C2O3
3325         'XO2           ', &    !< XO2
3326         'XO2N          ', &    !< XO2N
3327         'cro           ', &    !< CRO
3328         'HNO3          ', &    !< HNO3
3329         'H2O2          ', &    !< H2O2
3330         'iso           ', &    !< ISO
3331         'ispd          ', &    !< ISPD
3332         'to2           ', &    !< TO2
3333         'open          ', &    !< OPEN
3334         'terp          ', &    !< TERP
3335         'ec_f          ', &    !< EC_f
3336         'ec_c          ', &    !< EC_c
3337         'pom_f         ', &    !< POM_f
3338         'pom_c         ', &    !< POM_c
3339         'ppm_f         ', &    !< PPM_f
3340         'ppm_c         ', &    !< PPM_c
3341         'na_ff         ', &    !< Na_ff
3342         'na_f          ', &    !< Na_f
3343         'na_c          ', &    !< Na_c
3344         'na_cc         ', &    !< Na_cc
3345         'na_ccc        ', &    !< Na_ccc
3346         'dust_ff       ', &    !< dust_ff
3347         'dust_f        ', &    !< dust_f
3348         'dust_c        ', &    !< dust_c
3349         'dust_cc       ', &    !< dust_cc
3350         'dust_ccc      ', &    !< dust_ccc
3351         'tpm10         ', &    !< tpm10
3352         'tpm25         ', &    !< tpm25
3353         'tss           ', &    !< tss
3354         'tdust         ', &    !< tdust
3355         'tc            ', &    !< tc
3356         'tcg           ', &    !< tcg
3357         'tsoa          ', &    !< tsoa
3358         'tnmvoc        ', &    !< tnmvoc
3359         'SOxa          ', &    !< SOxa
3360         'NOya          ', &    !< NOya
3361         'NHxa          ', &    !< NHxa
3362         'NO2_obs       ', &    !< NO2_obs
3363         'tpm10_biascorr', &    !< tpm10_biascorr
3364         'tpm25_biascorr', &    !< tpm25_biascorr
3365         'O3_biascorr   ' /)    !< o3_biascorr
3366!
3367!-- tracer mole mass:
3368    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
3369         xm_O * 2 + xm_N, &                         !< NO2
3370         xm_O + xm_N, &                             !< NO
3371         xm_O * 3, &                                !< O3
3372         xm_C + xm_O, &                             !< CO
3373         xm_H * 2 + xm_C + xm_O, &                  !< FORM
3374         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
3375         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
3376         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
3377         xm_H * 3 + xm_C, &                         !< PAR
3378         xm_H * 3 + xm_C * 2, &                     !< OLE
3379         xm_H * 4 + xm_C * 2, &                     !< ETH
3380         xm_H * 8 + xm_C * 7, &                     !< TOL
3381         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
3382         xm_H * 10 + xm_C * 8, &                    !< XYL
3383         xm_S + xm_O * 4, &                         !< SO4a_f
3384         xm_S + xm_O * 2, &                         !< SO2
3385         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
3386         xm_H * 4 + xm_C, &                         !< CH4
3387         xm_H * 3 + xm_N, &                         !< NH3
3388         xm_O * 3 + xm_N, &                         !< NO3
3389         xm_H + xm_O, &                             !< OH
3390         xm_H + xm_O * 2, &                         !< HO2
3391         xm_O * 5 + xm_N * 2, &                     !< N2O5
3392         xm_S + xm_O * 4, &                         !< SO4a_c
3393         xm_H * 4 + xm_N, &                         !< NH4a_f
3394         xm_O * 3 + xm_N, &                         !< NO3a_f
3395         xm_O * 3 + xm_N, &                         !< NO3a_c
3396         xm_C * 2 + xm_O * 3, &                     !< C2O3
3397         xm_dummy, &                                !< XO2
3398         xm_dummy, &                                !< XO2N
3399         xm_dummy, &                                !< CRO
3400         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
3401         xm_H * 2 + xm_O * 2, &                     !< H2O2
3402         xm_H * 8 + xm_C * 5, &                     !< ISO
3403         xm_dummy, &                                !< ISPD
3404         xm_dummy, &                                !< TO2
3405         xm_dummy, &                                !< OPEN
3406         xm_H * 16 + xm_C * 10, &                   !< TERP
3407         xm_dummy, &                                !< EC_f
3408         xm_dummy, &                                !< EC_c
3409         xm_dummy, &                                !< POM_f
3410         xm_dummy, &                                !< POM_c
3411         xm_dummy, &                                !< PPM_f
3412         xm_dummy, &                                !< PPM_c
3413         xm_Na, &                                   !< Na_ff
3414         xm_Na, &                                   !< Na_f
3415         xm_Na, &                                   !< Na_c
3416         xm_Na, &                                   !< Na_cc
3417         xm_Na, &                                   !< Na_ccc
3418         xm_dummy, &                                !< dust_ff
3419         xm_dummy, &                                !< dust_f
3420         xm_dummy, &                                !< dust_c
3421         xm_dummy, &                                !< dust_cc
3422         xm_dummy, &                                !< dust_ccc
3423         xm_dummy, &                                !< tpm10
3424         xm_dummy, &                                !< tpm25
3425         xm_dummy, &                                !< tss
3426         xm_dummy, &                                !< tdust
3427         xm_dummy, &                                !< tc
3428         xm_dummy, &                                !< tcg
3429         xm_dummy, &                                !< tsoa
3430         xm_dummy, &                                !< tnmvoc
3431         xm_dummy, &                                !< SOxa
3432         xm_dummy, &                                !< NOya
3433         xm_dummy, &                                !< NHxa
3434         xm_O * 2 + xm_N, &                         !< NO2_obs
3435         xm_dummy, &                                !< tpm10_biascorr
3436         xm_dummy, &                                !< tpm25_biascorr
3437         xm_O * 3 /)                                !< o3_biascorr
3438!
3439!-- Initialize surface element m
3440    m = 0
3441    k = 0
3442!
3443!-- LSM or USM surface present at i,j:
3444!-- Default surfaces are NOT considered for deposition
3445    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
3446    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
3447!
3448!--For LSM surfaces
3449
3450    IF ( match_lsm )  THEN
3451!
3452!--    Get surface element information at i,j:
3453       m = surf_lsm_h%start_index(j,i)
3454       k = surf_lsm_h%k(m)
3455!
3456!--    Get needed variables for surface element m
3457       ustar_surf  = surf_lsm_h%us(m)
3458       z0h_surf    = surf_lsm_h%z0h(m)
3459       r_aero_surf = surf_lsm_h%r_a(m)
3460       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
3461       
3462       lai = surf_lsm_h%lai(m)
3463       sai = lai + 1
3464!
3465!--    For small grid spacing neglect R_a
3466       IF ( dzw(k) <= 1.0 )  THEN
3467          r_aero_surf = 0.0_wp
3468       ENDIF
3469!
3470!--    Initialize lu's
3471       luv_palm = 0
3472       luv_dep = 0
3473       lup_palm = 0
3474       lup_dep = 0
3475       luw_palm = 0
3476       luw_dep = 0
3477!
3478!--    Initialize budgets
3479       bud_luv = 0.0_wp
3480       bud_lup = 0.0_wp
3481       bud_luw = 0.0_wp
3482!
3483!--    Get land use for i,j and assign to DEPAC lu
3484       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3485          luv_palm = surf_lsm_h%vegetation_type(m)
3486          IF ( luv_palm == ind_luv_user )  THEN
3487             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3488             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
3489          ELSEIF ( luv_palm == ind_luv_b_soil )  THEN
3490             luv_dep = 9
3491          ELSEIF ( luv_palm == ind_luv_mixed_crops )  THEN
3492             luv_dep = 2
3493          ELSEIF ( luv_palm == ind_luv_s_grass )  THEN
3494             luv_dep = 1
3495          ELSEIF ( luv_palm == ind_luv_ev_needle_trees )  THEN
3496             luv_dep = 4
3497          ELSEIF ( luv_palm == ind_luv_de_needle_trees )  THEN
3498             luv_dep = 4
3499          ELSEIF ( luv_palm == ind_luv_ev_broad_trees )  THEN
3500             luv_dep = 12
3501          ELSEIF ( luv_palm == ind_luv_de_broad_trees )  THEN
3502             luv_dep = 5
3503          ELSEIF ( luv_palm == ind_luv_t_grass )  THEN
3504             luv_dep = 1
3505          ELSEIF ( luv_palm == ind_luv_desert )  THEN
3506             luv_dep = 9
3507          ELSEIF ( luv_palm == ind_luv_tundra )  THEN
3508             luv_dep = 8
3509          ELSEIF ( luv_palm == ind_luv_irr_crops )  THEN
3510             luv_dep = 2
3511          ELSEIF ( luv_palm == ind_luv_semidesert )  THEN
3512             luv_dep = 8
3513          ELSEIF ( luv_palm == ind_luv_ice )  THEN
3514             luv_dep = 10
3515          ELSEIF ( luv_palm == ind_luv_marsh )  THEN
3516             luv_dep = 8
3517          ELSEIF ( luv_palm == ind_luv_ev_shrubs )  THEN
3518             luv_dep = 14
3519          ELSEIF ( luv_palm == ind_luv_de_shrubs )  THEN
3520             luv_dep = 14
3521          ELSEIF ( luv_palm == ind_luv_mixed_forest )  THEN
3522             luv_dep = 4
3523          ELSEIF ( luv_palm == ind_luv_intrup_forest )  THEN
3524             luv_dep = 8     
3525          ENDIF
3526       ENDIF
3527
3528       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3529          lup_palm = surf_lsm_h%pavement_type(m)
3530          IF ( lup_palm == ind_lup_user )  THEN
3531             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3532             CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
3533          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
3534             lup_dep = 9
3535          ELSEIF ( lup_palm == ind_lup_asph )  THEN
3536             lup_dep = 9
3537          ELSEIF ( lup_palm == ind_lup_conc )  THEN
3538             lup_dep = 9
3539          ELSEIF ( lup_palm == ind_lup_sett )  THEN
3540             lup_dep = 9
3541          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
3542             lup_dep = 9
3543          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
3544             lup_dep = 9       
3545          ELSEIF ( lup_palm == ind_lup_metal )  THEN
3546             lup_dep = 9
3547          ELSEIF ( lup_palm == ind_lup_wood )  THEN
3548             lup_dep = 9   
3549          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
3550             lup_dep = 9
3551          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
3552             lup_dep = 9
3553          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
3554             lup_dep = 9
3555          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
3556             lup_dep = 9
3557          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
3558             lup_dep = 9
3559          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
3560             lup_dep = 9
3561          ELSEIF ( lup_palm == ind_lup_clay )  THEN
3562             lup_dep = 9
3563          ENDIF
3564       ENDIF
3565
3566       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3567          luw_palm = surf_lsm_h%water_type(m)     
3568          IF ( luw_palm == ind_luw_user )  THEN
3569             message_string = 'No water type defined. Please define water type to enable deposition calculation'
3570             CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
3571          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
3572             luw_dep = 13
3573          ELSEIF ( luw_palm == ind_luw_river )  THEN
3574             luw_dep = 13
3575          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
3576             luw_dep = 6
3577          ELSEIF ( luw_palm == ind_luw_pond )  THEN
3578             luw_dep = 13
3579          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
3580             luw_dep = 13 
3581          ENDIF
3582       ENDIF
3583!
3584!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
3585       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
3586          nwet = 1
3587       ELSE
3588          nwet = 0
3589       ENDIF
3590!
3591!--    Compute length of time step
3592       IF ( call_chem_at_all_substeps )  THEN
3593          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3594       ELSE
3595          dt_chem = dt_3d
3596       ENDIF
3597
3598       dh = dzw(k)
3599       inv_dh = 1.0_wp / dh
3600       dt_dh = dt_chem/dh
3601!
3602!--    Concentration at i,j,k
3603       DO  lsp = 1, nspec
3604          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3605       ENDDO
3606
3607!--    Temperature at i,j,k
3608       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3609
3610       ts       = temp_tmp - 273.15  !< in degrees celcius
3611!
3612!--    Viscosity of air
3613       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
3614!
3615!--    Air density at k
3616       dens = rho_air_zw(k)
3617!
3618!--    Calculate relative humidity from specific humidity for DEPAC
3619       qv_tmp = MAX(q(k,j,i),0.0_wp)
3620       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
3621!
3622!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3623!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3624!
3625!--    Vegetation
3626       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3627
3628!
3629!--       No vegetation on bare soil, desert or ice:
3630          IF ( ( luv_palm == ind_luv_b_soil ) .OR. &
3631                ( luv_palm == ind_luv_desert ) .OR. &
3632                 ( luv_palm == ind_luv_ice ) ) THEN
3633
3634             lai = 0.0_wp
3635             sai = 0.0_wp
3636
3637          ENDIF
3638         
3639          slinnfac = 1.0_wp
3640!
3641!--       Get deposition velocity vd
3642          DO  lsp = 1, nvar
3643!
3644!--          Initialize
3645             vs = 0.0_wp
3646             vd_lu = 0.0_wp
3647             rs = 0.0_wp
3648             rb = 0.0_wp
3649             rc_tot = 0.0_wp
3650             IF ( spc_names(lsp) == 'PM10' )  THEN
3651                part_type = 1
3652!
3653!--             Sedimentation velocity
3654                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3655                     particle_pars(ind_p_size, part_type), &
3656                     particle_pars(ind_p_slip, part_type), &
3657                     visc)
3658
3659                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3660                     vs, &
3661                     particle_pars(ind_p_size, part_type), &
3662                     particle_pars(ind_p_slip, part_type), &
3663                     nwet, temp_tmp, dens, visc, &
3664                     luv_dep,  &
3665                     r_aero_surf, ustar_surf )
3666
3667                bud_luv(lsp) = - conc_ijk(lsp) * &
3668                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3669
3670
3671             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3672                part_type = 2
3673!
3674!--             Sedimentation velocity
3675                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3676                     particle_pars(ind_p_size, part_type), &
3677                     particle_pars(ind_p_slip, part_type), &
3678                     visc)
3679
3680                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3681                     vs, &
3682                     particle_pars(ind_p_size, part_type), &
3683                     particle_pars(ind_p_slip, part_type), &
3684                     nwet, temp_tmp, dens, visc, &
3685                     luv_dep , &
3686                     r_aero_surf, ustar_surf )
3687
3688                bud_luv(lsp) = - conc_ijk(lsp) * &
3689                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3690
3691             ELSE !< GASES
3692!
3693!--             Read spc_name of current species for gas parameter
3694                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
3695                   i_pspec = 0
3696                   DO  pspec = 1, nposp
3697                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3698                         i_pspec = pspec
3699                      END IF
3700                   ENDDO
3701
3702                ELSE
3703!
3704!--             For now species not deposited
3705                   CYCLE
3706                ENDIF
3707!
3708!--             Factor used for conversion from ppb to ug/m3 :
3709!--             ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3710!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3711!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
3712!--             thus:
3713!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3714!--             Use density at k:
3715
3716                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3717!
3718!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3719                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3720                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3721!
3722!--             Diffusivity for DEPAC relevant gases
3723!--             Use default value
3724                diffusivity            = 0.11e-4
3725!
3726!--             overwrite with known coefficients of diffusivity from Massman (1998)
3727                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3728                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3729                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3730                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3731                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3732                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3733                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3734!
3735!--             Get quasi-laminar boundary layer resistance rb:
3736                CALL get_rb_cell( (luv_dep == ilu_water_sea) .OR. (luv_dep == ilu_water_inland), &
3737                     z0h_surf, ustar_surf, diffusivity, &
3738                     rb )
3739!
3740!--             Get rc_tot
3741                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
3742                     rh_surf, lai, sai, nwet, luv_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3743                     r_aero_surf , rb )
3744!
3745!--             Calculate budget
3746                IF ( rc_tot <= 0.0 )  THEN
3747
3748                   bud_luv(lsp) = 0.0_wp
3749
3750                ELSE
3751
3752                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3753
3754                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3755                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3756                ENDIF
3757
3758             ENDIF
3759          ENDDO
3760       ENDIF
3761!
3762!--    Pavement
3763       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3764!
3765!--       No vegetation on pavements:
3766          lai = 0.0_wp
3767          sai = 0.0_wp
3768
3769          slinnfac = 1.0_wp
3770!
3771!--       Get vd
3772          DO  lsp = 1, nvar
3773!
3774!--       Initialize
3775             vs = 0.0_wp
3776             vd_lu = 0.0_wp
3777             rs = 0.0_wp
3778             rb = 0.0_wp
3779             rc_tot = 0.0_wp
3780             IF ( spc_names(lsp) == 'PM10' )  THEN
3781                part_type = 1
3782!
3783!--             Sedimentation velocity:
3784                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3785                     particle_pars(ind_p_size, part_type), &
3786                     particle_pars(ind_p_slip, part_type), &
3787                     visc)
3788
3789                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3790                     vs, &
3791                     particle_pars(ind_p_size, part_type), &
3792                     particle_pars(ind_p_slip, part_type), &
3793                     nwet, temp_tmp, dens, visc, &
3794                     lup_dep,  &
3795                     r_aero_surf, ustar_surf )
3796
3797                bud_lup(lsp) = - conc_ijk(lsp) * &
3798                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3799
3800
3801             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3802                part_type = 2
3803!
3804!--             Sedimentation velocity:
3805                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3806                     particle_pars(ind_p_size, part_type), &
3807                     particle_pars(ind_p_slip, part_type), &
3808                     visc)
3809
3810                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3811                     vs, &
3812                     particle_pars(ind_p_size, part_type), &
3813                     particle_pars(ind_p_slip, part_type), &
3814                     nwet, temp_tmp, dens, visc, &
3815                     lup_dep, &
3816                     r_aero_surf, ustar_surf )
3817
3818                bud_lup(lsp) = - conc_ijk(lsp) * &
3819                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3820
3821             ELSE  !<GASES
3822!
3823!--             Read spc_name of current species for gas parameter
3824
3825                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3826                   i_pspec = 0
3827                   DO  pspec = 1, nposp
3828                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3829                         i_pspec = pspec
3830                      END IF
3831                   ENDDO
3832
3833                ELSE
3834!
3835!--                For now species not deposited
3836                   CYCLE
3837                ENDIF
3838!
3839!--             Factor used for conversion from ppb to ug/m3 :
3840!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3841!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3842!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3843!--             thus:
3844!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3845!--             Use density at lowest layer:
3846
3847                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3848!
3849!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3850                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3851                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3852!
3853!--             Diffusivity for DEPAC relevant gases
3854!--             Use default value
3855                diffusivity            = 0.11e-4
3856!
3857!--             overwrite with known coefficients of diffusivity from Massman (1998)
3858                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3859                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3860                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3861                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3862                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3863                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3864                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3865!
3866!--             Get quasi-laminar boundary layer resistance rb:
3867                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
3868                     z0h_surf, ustar_surf, diffusivity, rb )
3869!
3870!--             Get rc_tot
3871                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
3872                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
3873                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
3874                                         r_aero_surf , rb )
3875!
3876!--             Calculate budget
3877                IF ( rc_tot <= 0.0 )  THEN
3878                   bud_lup(lsp) = 0.0_wp
3879                ELSE
3880                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3881                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3882                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3883                ENDIF
3884
3885             ENDIF
3886          ENDDO
3887       ENDIF
3888!
3889!--    Water
3890       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3891!
3892!--       No vegetation on water:
3893          lai = 0.0_wp
3894          sai = 0.0_wp
3895          slinnfac = 1.0_wp
3896!
3897!--       Get vd
3898          DO  lsp = 1, nvar
3899!
3900!--          Initialize
3901             vs = 0.0_wp
3902             vd_lu = 0.0_wp
3903             rs = 0.0_wp
3904             rb = 0.0_wp
3905             rc_tot = 0.0_wp 
3906             IF ( spc_names(lsp) == 'PM10' )  THEN
3907                part_type = 1
3908!
3909!--             Sedimentation velocity:
3910                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3911                     particle_pars(ind_p_size, part_type), &
3912                     particle_pars(ind_p_slip, part_type), &
3913                     visc)
3914
3915                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3916                     vs, &
3917                     particle_pars(ind_p_size, part_type), &
3918                     particle_pars(ind_p_slip, part_type), &
3919                     nwet, temp_tmp, dens, visc, &
3920                     luw_dep, &
3921                     r_aero_surf, ustar_surf )
3922
3923                bud_luw(lsp) = - conc_ijk(lsp) * &
3924                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3925
3926             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3927                part_type = 2
3928!
3929!--             Sedimentation velocity:
3930                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3931                     particle_pars(ind_p_size, part_type), &
3932                     particle_pars(ind_p_slip, part_type), &
3933                     visc)
3934
3935                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3936                     vs, &
3937                     particle_pars(ind_p_size, part_type), &
3938                     particle_pars(ind_p_slip, part_type), &
3939                     nwet, temp_tmp, dens, visc, &
3940                     luw_dep, &
3941                     r_aero_surf, ustar_surf )
3942
3943                bud_luw(lsp) = - conc_ijk(lsp) * &
3944                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3945
3946             ELSE  !<GASES
3947!
3948!--             Read spc_name of current species for gas PARAMETER
3949
3950                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3951                   i_pspec = 0
3952                   DO  pspec = 1, nposp
3953                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3954                         i_pspec = pspec
3955                      END IF
3956                   ENDDO
3957
3958                ELSE
3959!
3960!--                For now species not deposited
3961                   CYCLE
3962                ENDIF
3963!
3964!--             Factor used for conversion from ppb to ug/m3 :
3965!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3966!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3967!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3968!--             thus:
3969!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3970!--             Use density at lowest layer:
3971
3972                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3973!
3974!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3975!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3976                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
3977!
3978!--             Diffusivity for DEPAC relevant gases
3979!--             Use default value
3980                diffusivity            = 0.11e-4
3981!
3982!--             overwrite with known coefficients of diffusivity from Massman (1998)
3983                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3984                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3985                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3986                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3987                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3988                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3989                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3990!
3991!--             Get quasi-laminar boundary layer resistance rb:
3992                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
3993                     z0h_surf, ustar_surf, diffusivity, rb )
3994
3995!--             Get rc_tot
3996                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         & 
3997                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
3998                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
3999                                         r_aero_surf , rb )
4000!
4001!--             Calculate budget
4002                IF ( rc_tot <= 0.0 )  THEN
4003
4004                   bud_luw(lsp) = 0.0_wp
4005
4006                ELSE
4007
4008                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4009
4010                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4011                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4012                ENDIF
4013
4014             ENDIF
4015          ENDDO
4016       ENDIF
4017
4018
4019       bud = 0.0_wp
4020!
4021!--    Calculate overall budget for surface m and adapt concentration
4022       DO  lsp = 1, nspec
4023
4024          bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
4025               surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
4026               surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
4027!
4028!--       Compute new concentration:
4029          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4030
4031          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
4032
4033       ENDDO
4034
4035    ENDIF
4036!
4037!-- For USM surfaces   
4038
4039    IF ( match_usm )  THEN
4040!
4041!--    Get surface element information at i,j:
4042       m = surf_usm_h%start_index(j,i)
4043       k = surf_usm_h%k(m)
4044!
4045!--    Get needed variables for surface element m
4046       ustar_surf  = surf_usm_h%us(m)
4047       z0h_surf    = surf_usm_h%z0h(m)
4048       r_aero_surf = surf_usm_h%r_a(m)
4049       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
4050       lai = surf_usm_h%lai(m)
4051       sai = lai + 1
4052!
4053!--    For small grid spacing neglect R_a
4054       IF ( dzw(k) <= 1.0 )  THEN
4055          r_aero_surf = 0.0_wp
4056       ENDIF
4057!
4058!--    Initialize lu's
4059       luu_palm = 0
4060       luu_dep = 0
4061       lug_palm = 0
4062       lug_dep = 0
4063       lud_palm = 0
4064       lud_dep = 0
4065!
4066!--    Initialize budgets
4067       bud_luu  = 0.0_wp
4068       bud_lug = 0.0_wp
4069       bud_lud = 0.0_wp
4070!
4071!--    Get land use for i,j and assign to DEPAC lu
4072       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4073!
4074!--       For green urban surfaces (e.g. green roofs
4075!--       assume LU short grass
4076          lug_palm = ind_luv_s_grass
4077          IF ( lug_palm == ind_luv_user )  THEN
4078             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
4079             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
4080          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
4081             lug_dep = 9
4082          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
4083             lug_dep = 2
4084          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
4085             lug_dep = 1
4086          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
4087             lug_dep = 4
4088          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
4089             lug_dep = 4
4090          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
4091             lug_dep = 12
4092          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
4093             lug_dep = 5
4094          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
4095             lug_dep = 1
4096          ELSEIF ( lug_palm == ind_luv_desert )  THEN
4097             lug_dep = 9
4098          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
4099             lug_dep = 8
4100          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
4101             lug_dep = 2
4102          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
4103             lug_dep = 8
4104          ELSEIF ( lug_palm == ind_luv_ice )  THEN
4105             lug_dep = 10
4106          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
4107             lug_dep = 8
4108          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
4109             lug_dep = 14
4110          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
4111             lug_dep = 14
4112          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
4113             lug_dep = 4
4114          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
4115             lug_dep = 8     
4116          ENDIF
4117       ENDIF
4118
4119       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
4120!
4121!--       For walls in USM assume concrete walls/roofs,
4122!--       assumed LU class desert as also assumed for
4123!--       pavements in LSM
4124          luu_palm = ind_lup_conc
4125          IF ( luu_palm == ind_lup_user )  THEN
4126             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
4127             CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
4128          ELSEIF ( luu_palm == ind_lup_asph_conc )  THEN
4129             luu_dep = 9
4130          ELSEIF ( luu_palm == ind_lup_asph )  THEN
4131             luu_dep = 9
4132          ELSEIF ( luu_palm ==  ind_lup_conc )  THEN
4133             luu_dep = 9
4134          ELSEIF ( luu_palm ==  ind_lup_sett )  THEN
4135             luu_dep = 9
4136          ELSEIF ( luu_palm == ind_lup_pav_stones )  THEN
4137             luu_dep = 9
4138          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
4139             luu_dep = 9       
4140          ELSEIF ( luu_palm == ind_lup_metal )  THEN
4141             luu_dep = 9
4142          ELSEIF ( luu_palm == ind_lup_wood )  THEN
4143             luu_dep = 9   
4144          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
4145             luu_dep = 9
4146          ELSEIF ( luu_palm == ind_lup_f_gravel )  THEN
4147             luu_dep = 9
4148          ELSEIF ( luu_palm == ind_lup_pebblest )  THEN
4149             luu_dep = 9
4150          ELSEIF ( luu_palm == ind_lup_woodchips )  THEN
4151             luu_dep = 9
4152          ELSEIF ( luu_palm == ind_lup_tartan )  THEN
4153             luu_dep = 9
4154          ELSEIF ( luu_palm == ind_lup_art_turf )  THEN
4155             luu_dep = 9
4156          ELSEIF ( luu_palm == ind_lup_clay )  THEN
4157             luu_dep = 9
4158          ENDIF
4159       ENDIF
4160
4161       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4162!
4163!--       For windows in USM assume metal as this is
4164!--       as close as we get, assumed LU class desert
4165!--       as also assumed for pavements in LSM
4166          lud_palm = ind_lup_metal     
4167          IF ( lud_palm == ind_lup_user )  THEN
4168             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
4169             CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
4170          ELSEIF ( lud_palm == ind_lup_asph_conc )  THEN
4171             lud_dep = 9
4172          ELSEIF ( lud_palm == ind_lup_asph )  THEN
4173             lud_dep = 9
4174          ELSEIF ( lud_palm ==  ind_lup_conc )  THEN
4175             lud_dep = 9
4176          ELSEIF ( lud_palm ==  ind_lup_sett )  THEN
4177             lud_dep = 9
4178          ELSEIF ( lud_palm == ind_lup_pav_stones )  THEN
4179             lud_dep = 9
4180          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
4181             lud_dep = 9       
4182          ELSEIF ( lud_palm == ind_lup_metal )  THEN
4183             lud_dep = 9
4184          ELSEIF ( lud_palm == ind_lup_wood )  THEN
4185             lud_dep = 9   
4186          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
4187             lud_dep = 9
4188          ELSEIF ( lud_palm == ind_lup_f_gravel )  THEN
4189             lud_dep = 9
4190          ELSEIF ( lud_palm == ind_lup_pebblest )  THEN
4191             lud_dep = 9
4192          ELSEIF ( lud_palm == ind_lup_woodchips )  THEN
4193             lud_dep = 9
4194          ELSEIF ( lud_palm == ind_lup_tartan )  THEN
4195             lud_dep = 9
4196          ELSEIF ( lud_palm == ind_lup_art_turf )  THEN
4197             lud_dep = 9
4198          ELSEIF ( lud_palm == ind_lup_clay )  THEN
4199             lud_dep = 9
4200          ENDIF
4201       ENDIF
4202!
4203!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
4204!--    Set wetness indicator to dry or wet for usm vegetation or pavement
4205       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
4206       !   nwet = 1
4207       !ELSE
4208       nwet = 0
4209       !ENDIF
4210!
4211!--    Compute length of time step
4212       IF ( call_chem_at_all_substeps )  THEN
4213          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
4214       ELSE
4215          dt_chem = dt_3d
4216       ENDIF
4217
4218       dh = dzw(k)
4219       inv_dh = 1.0_wp / dh
4220       dt_dh = dt_chem/dh
4221!
4222!--    Concentration at i,j,k
4223       DO  lsp = 1, nspec
4224          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
4225       ENDDO
4226!
4227!--    Temperature at i,j,k
4228       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
4229
4230       ts       = temp_tmp - 273.15  !< in degrees celcius
4231!
4232!--    Viscosity of air
4233       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
4234!
4235!--    Air density at k
4236       dens = rho_air_zw(k)
4237!
4238!--    Calculate relative humidity from specific humidity for DEPAC
4239       qv_tmp = MAX(q(k,j,i),0.0_wp)
4240       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
4241!
4242!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
4243!--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
4244!
4245!--    Walls/roofs
4246       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
4247!
4248!--       No vegetation on non-green walls:
4249          lai = 0.0_wp
4250          sai = 0.0_wp
4251
4252          slinnfac = 1.0_wp
4253!
4254!--       Get vd
4255          DO  lsp = 1, nvar
4256!
4257!--          Initialize
4258             vs = 0.0_wp
4259             vd_lu = 0.0_wp
4260             rs = 0.0_wp
4261             rb = 0.0_wp
4262             rc_tot = 0.0_wp
4263             IF (spc_names(lsp) == 'PM10' )  THEN
4264                part_type = 1
4265!
4266!--             Sedimentation velocity
4267                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4268                     particle_pars(ind_p_size, part_type), &
4269                     particle_pars(ind_p_slip, part_type), &
4270                     visc)
4271
4272                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4273                     vs, &
4274                     particle_pars(ind_p_size, part_type), &
4275                     particle_pars(ind_p_slip, part_type), &
4276                     nwet, temp_tmp, dens, visc, &
4277                     luu_dep,  &
4278                     r_aero_surf, ustar_surf )
4279
4280                bud_luu(lsp) = - conc_ijk(lsp) * &
4281                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4282
4283             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4284                part_type = 2
4285!
4286!--             Sedimentation velocity
4287                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4288                     particle_pars(ind_p_size, part_type), &
4289                     particle_pars(ind_p_slip, part_type), &
4290                     visc)
4291
4292                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4293                     vs, &
4294                     particle_pars(ind_p_size, part_type), &
4295                     particle_pars(ind_p_slip, part_type), &
4296                     nwet, temp_tmp, dens, visc, &
4297                     luu_dep , &
4298                     r_aero_surf, ustar_surf )
4299
4300                bud_luu(lsp) = - conc_ijk(lsp) * &
4301                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4302
4303             ELSE  !< GASES
4304!
4305!--             Read spc_name of current species for gas parameter
4306
4307                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4308                   i_pspec = 0
4309                   DO  pspec = 1, nposp
4310                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4311                         i_pspec = pspec
4312                      END IF
4313                   ENDDO
4314                ELSE
4315!
4316!--                For now species not deposited
4317                   CYCLE
4318                ENDIF
4319!
4320!--             Factor used for conversion from ppb to ug/m3 :
4321!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4322!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4323!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
4324!--             thus:
4325!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4326!--             Use density at k:
4327
4328                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4329
4330                !
4331!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4332!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4333                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4334!
4335!--             Diffusivity for DEPAC relevant gases
4336!--             Use default value
4337                diffusivity            = 0.11e-4
4338!
4339!--             overwrite with known coefficients of diffusivity from Massman (1998)
4340                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4341                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4342                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4343                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4344                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4345                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4346                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4347!
4348!--             Get quasi-laminar boundary layer resistance rb:
4349                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
4350                     z0h_surf, ustar_surf, diffusivity, &
4351                     rb )
4352!
4353!--             Get rc_tot
4354                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
4355                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
4356                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
4357                                         r_aero_surf, rb )
4358!
4359!--             Calculate budget
4360                IF ( rc_tot <= 0.0 )  THEN
4361
4362                   bud_luu(lsp) = 0.0_wp
4363
4364                ELSE
4365
4366                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4367
4368                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4369                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4370                ENDIF
4371
4372             ENDIF
4373          ENDDO
4374       ENDIF
4375!
4376!--    Green usm surfaces
4377       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4378
4379!
4380!--       No vegetation on bare soil, desert or ice:
4381          IF ( ( lug_palm == ind_luv_b_soil ) .OR. &
4382                ( lug_palm == ind_luv_desert ) .OR. &
4383                 ( lug_palm == ind_luv_ice ) ) THEN
4384
4385             lai = 0.0_wp
4386             sai = 0.0_wp
4387
4388          ENDIF
4389
4390         
4391          slinnfac = 1.0_wp
4392!
4393!--       Get vd
4394          DO  lsp = 1, nvar
4395!
4396!--          Initialize
4397             vs = 0.0_wp
4398             vd_lu = 0.0_wp
4399             rs = 0.0_wp
4400             rb = 0.0_wp
4401             rc_tot = 0.0_wp
4402             IF ( spc_names(lsp) == 'PM10' )  THEN
4403                part_type = 1
4404!
4405!--             Sedimentation velocity
4406                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4407                     particle_pars(ind_p_size, part_type), &
4408                     particle_pars(ind_p_slip, part_type), &
4409                     visc)
4410
4411                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4412                     vs, &
4413                     particle_pars(ind_p_size, part_type), &
4414                     particle_pars(ind_p_slip, part_type), &
4415                     nwet, temp_tmp, dens, visc, &
4416                     lug_dep,  &
4417                     r_aero_surf, ustar_surf )
4418
4419                bud_lug(lsp) = - conc_ijk(lsp) * &
4420                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4421
4422             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4423                part_type = 2
4424!
4425!--             Sedimentation velocity
4426                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4427                     particle_pars(ind_p_size, part_type), &
4428                     particle_pars(ind_p_slip, part_type), &
4429                     visc)
4430
4431                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4432                     vs, &
4433                     particle_pars(ind_p_size, part_type), &
4434                     particle_pars(ind_p_slip, part_type), &
4435                     nwet, temp_tmp, dens, visc, &
4436                     lug_dep, &
4437                     r_aero_surf, ustar_surf )
4438
4439                bud_lug(lsp) = - conc_ijk(lsp) * &
4440                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4441
4442             ELSE  !< GASES
4443!
4444!--             Read spc_name of current species for gas parameter
4445
4446                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4447                   i_pspec = 0
4448                   DO  pspec = 1, nposp
4449                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4450                         i_pspec = pspec
4451                      END IF
4452                   ENDDO
4453                ELSE
4454!
4455!--                For now species not deposited
4456                   CYCLE
4457                ENDIF
4458!
4459!--             Factor used for conversion from ppb to ug/m3 :
4460!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4461!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4462!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4463!--             thus:
4464!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4465!--             Use density at k:
4466
4467                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4468!
4469!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4470                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4471                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4472!
4473!--             Diffusivity for DEPAC relevant gases
4474!--             Use default value
4475                diffusivity            = 0.11e-4
4476!
4477!--             overwrite with known coefficients of diffusivity from Massman (1998)
4478                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4479                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4480                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4481                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4482                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4483                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4484                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4485!
4486!--             Get quasi-laminar boundary layer resistance rb:
4487                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
4488                     z0h_surf, ustar_surf, diffusivity, &
4489                     rb )
4490!
4491!--             Get rc_tot
4492                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
4493                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
4494                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
4495                                         r_aero_surf , rb )
4496!
4497!--             Calculate budget
4498                IF ( rc_tot <= 0.0 )  THEN
4499
4500                   bud_lug(lsp) = 0.0_wp
4501
4502                ELSE
4503
4504                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4505
4506                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4507                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4508                ENDIF
4509
4510             ENDIF
4511          ENDDO
4512       ENDIF
4513!
4514!--    Windows
4515       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4516!
4517!--       No vegetation on windows:
4518          lai = 0.0_wp
4519          sai = 0.0_wp
4520
4521          slinnfac = 1.0_wp
4522!
4523!--       Get vd
4524          DO  lsp = 1, nvar
4525!
4526!--          Initialize
4527             vs = 0.0_wp
4528             vd_lu = 0.0_wp
4529             rs = 0.0_wp
4530             rb = 0.0_wp
4531             rc_tot = 0.0_wp 
4532             IF ( spc_names(lsp) == 'PM10' )  THEN
4533                part_type = 1
4534!
4535!--             Sedimentation velocity
4536                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4537                     particle_pars(ind_p_size, part_type), &
4538                     particle_pars(ind_p_slip, part_type), &
4539                     visc)
4540
4541                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4542                     particle_pars(ind_p_size, part_type), &
4543                     particle_pars(ind_p_slip, part_type), &
4544                     nwet, temp_tmp, dens, visc,              &
4545                     lud_dep, r_aero_surf, ustar_surf )
4546
4547                bud_lud(lsp) = - conc_ijk(lsp) * &
4548                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4549
4550             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4551                part_type = 2
4552!
4553!--             Sedimentation velocity
4554                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4555                     particle_pars(ind_p_size, part_type), &
4556                     particle_pars(ind_p_slip, part_type), &
4557                     visc)
4558
4559                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4560                     particle_pars(ind_p_size, part_type), &
4561                     particle_pars(ind_p_slip, part_type), &
4562                     nwet, temp_tmp, dens, visc, &
4563                     lud_dep, &
4564                     r_aero_surf, ustar_surf )
4565
4566                bud_lud(lsp) = - conc_ijk(lsp) * &
4567                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4568
4569             ELSE  !< GASES
4570!
4571!--             Read spc_name of current species for gas PARAMETER
4572
4573                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4574                   i_pspec = 0
4575                   DO  pspec = 1, nposp
4576                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4577                         i_pspec = pspec
4578                      END IF
4579                   ENDDO
4580                ELSE
4581!
4582!--                For now species not deposited
4583                   CYCLE
4584                ENDIF
4585!
4586!--             Factor used for conversion from ppb to ug/m3 :
4587!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4588!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4589!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4590!--             thus:
4591!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4592!--             Use density at k:
4593
4594                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4595!
4596!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4597!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4598                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4599!
4600!--             Diffusivity for DEPAC relevant gases
4601!--             Use default value
4602                diffusivity = 0.11e-4
4603!
4604!--             overwrite with known coefficients of diffusivity from Massman (1998)
4605                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4606                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4607                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4608                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4609                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4610                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4611                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4612!
4613!--             Get quasi-laminar boundary layer resistance rb:
4614                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
4615                     z0h_surf, ustar_surf, diffusivity, rb )
4616!
4617!--             Get rc_tot
4618                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
4619                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
4620                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
4621                                         r_aero_surf , rb )
4622!
4623!--             Calculate budget
4624                IF ( rc_tot <= 0.0 )  THEN
4625
4626                   bud_lud(lsp) = 0.0_wp
4627
4628                ELSE
4629
4630                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4631
4632                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4633                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4634                ENDIF
4635
4636             ENDIF
4637          ENDDO
4638       ENDIF
4639
4640
4641       bud = 0.0_wp
4642!
4643!--    Calculate overall budget for surface m and adapt concentration
4644       DO  lsp = 1, nspec
4645
4646
4647          bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
4648               surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
4649               surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
4650!
4651!--       Compute new concentration
4652          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4653
4654          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
4655
4656       ENDDO
4657
4658    ENDIF
4659
4660
4661 END SUBROUTINE chem_depo
4662
4663
4664!------------------------------------------------------------------------------!
4665! Description:
4666! ------------
4667!> Subroutine to compute total canopy (or surface) resistance Rc for gases
4668!>
4669!> DEPAC:
4670!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
4671!> module of the LOTOS-EUROS model (Manders et al., 2017)
4672!>
4673!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
4674!> van Zanten et al., 2010.
4675!------------------------------------------------------------------------------!
4676 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
4677      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
4678      ra, rb ) 
4679!
4680!--   Some of depac arguments are OPTIONAL:
4681!--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
4682!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
4683!--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
4684!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4685!--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
4686!--
4687!--    C. compute effective Rc based on compensation points (used for OPS):
4688!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4689!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4690!--               ra, rb, rc_eff)
4691!--    X1. Extra (OPTIONAL) output variables:
4692!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4693!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4694!--               ra, rb, rc_eff, &