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

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

Bugfix in chemistry decycling; Replace global arrays to setup chemical emissions by local ones

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