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

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