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

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

Merge from branch resler: Changed behaviour of masked output over surface to follow terrain and ignore buildings

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