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

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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