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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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