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

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

Moved loop over chem_species into chem_boundary_conds_decycle

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