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

Last change on this file since 3876 was 3876, checked in by knoop, 2 years ago

Moved "photolysis_scheme", "chem_species" and "phot_frequen" to chem_modules

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