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

Last change on this file since 4290 was 4290, checked in by banzhafs, 17 months ago

Bugfix in sedimentation resistance calculation in drydepo_aero_zhang_vd subroutine

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