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

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

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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