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

Last change on this file since 4306 was 4306, checked in by banzhafs, 19 months ago

Correction to r4304 commit in chemistry_model_mod

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