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

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

bugfixes: uninitialized variable in dry deposition; emission output

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