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

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

Bugfix in setting flags inidicating wall-bounded grid-points (used for control TKE production near walls); Bugfix in chemistry decycling initialization

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